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Abstract 

In this paper, a new re-initialization method for the conservative level-set func¬ 
tion is put forward. First, it has been shown that the re-initialization and 
advection equations of the conservative level-set function are mathematically 
equivalent to the re-initialization and advection equations of the localized signed 
distance function. Next, a new discretization for the spatial derivatives of the 
conservative level-set function has been proposed. This new discretization is 
consistent with the re-initialization procedure and it guarantees a second-order 
convergence rate of the interface curvature on gradually refined grids. The new 
re-initialization method does not introduce artificial deformations to stationary 
and non-stationary interfaces, even when the number of re-initialization steps 
is large. 
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1. Introduction 

In the conservative level-set method introduced by Olsson and Kreiss [1], 
a solution of the transport equation of the characteristic function a (x, t) is 
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divided into two steps: advection and re-initialization. These two steps are 
carried out consecutively in only one single time step At. The solution of the 
a (x, t) advection equation is typically performed using an explicit or implicit 
time discretization with a TVD MUSCL flux limiter, in order to keep 0 < 
cc (x, t) < la bounded function [1, 2, 3, 4, 5]. The present work provides 
improvements in the second step of the above interface capturing procedure. 
Let us consider the re-initialization equation 

= V-[D\Va\n r -Ca(l-a)n r ], ( 1 ) 

where a (x, t) denotes a regularized Heaviside function, r is the artificial time, 
nr = Va/|Va| is a vector normal to the iso-lines (iso-surfaces) of a (x, t). We 
assume that C = 1 m/s, D = Ax/2C = €hC where Ax = Ay = Az = 1/N C is 
the size a control volume and N c is the total number of control volumes in the 
x, y or z direction. 

In [1] it has been shown the solution of Eq. (1) after advection of a(x, f), 
reduces artificial deformations of the interface induced by numerical errors dur¬ 
ing the advection step. Therein it was also noticed that the analytical solution 
to the stationary equation (1) is given by 


a(V>o(x,f)) = 1- 


1 

l+exp(^ 0 CM) /f-h) 


1 

2 


1 + tanh 


( j>o (*, t) \' 

V 2 e h )y 


( 2 ) 


where il> 0 (x,t) is the signed distance function [5, 6, 7]. When e h —> 0, then 
a {ip o) in equation (2) tends to the phase indicator function which is given by 
the exact Heaviside function H (ipo) as is demonstrated in Appendix A. The 
phase indicator function H {ipo) is typically discretized in the volume of fluid 
(VOF) family of methods satisfying the law of conservation of mass [8]. In 
real simulations eh ^ 0, and hence a {ipo) in equation (2) is not represented 
by a sharp jump localized at the interface T. The level-set function a ( ipo ) 
is a Lipschitz continuous function and, therefore, resembles the signed distance 
function in the standard level-set (LS) method introduced by Oslier and Sethian 
[9] and extended by Sussman et al. [10, 11]. For these reasons, Olsson and 
Kreiss [1] called their interface capturing technique the conservative level-set 
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(CLS) method. The signed distance function derived from Eq. (2) is given by 
the equation 

i’o (a) = e h In , ■ (3) 

[1 - a (V>o)J 

We note that in Eq. (2) the interface T is localized at a (xr, t) = 1/2 whereas in 
Eq. (3) the set of points where the signed distance function ipo (a (xp, t) — 1/2) = 
0 represents a position of the interface E. 

A mapping between a {ip o) and ipo (a) in equations (2) and (3) suggests a 
closer relation between the CLS method and the standard LS method exists. In 
the present paper, we show how this relation can be established. Moreover, we 
put forward a new method for computation of higher-order spatial derivatives 
of a {ipo), which is consistent with the new re-initialization procedure. The 
spatial derivatives of a {ipo) obtained with our new method are later used to 
approximate the interface curvature n with second-order accuracy. 

A relation between the regularized Heaviside function and the signed dis¬ 
tance function was first observed by Glasner [6], and was used for a non-linear 
preconditioning of the phase-field equation. The non-linear preconditioning was 
later exploited by Sun and Beckermann [7] in order to solve the phase-field 
equation in a context of the interface capturing. Therein, it was mentioned 
that the stationary solutions to the phase-field equation and to equation (1) are 
different. Later in this paper, the key differences between the present results 
and the results reported in [7] are addressed. 

The main difficulty in using a {ipo) and ipo (a) interchangeably is the lack 
of a correct numerical solution to the re-initialization equation (1). Although 
vast literature concerning the numerical solution of the re-initialization equa¬ 
tion for the signed distance function exists, see [12, 13, 14] to name only a few 
recent works, the solution of equation (1) has drawn less attention. In [1] the 
re-initialization equation (1) is solved directly, in [15] to reduce artificial inter¬ 
face deformations due to discretization errors; the diffusive term in Eq. (1) was 
projected on nr, however as it is shown by Shukla et al. [2], this reformulation 
leads to numerical instabilities. Recently, McCaslin and Desjardins [16] pro- 
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posed to multiply diffusive and compressive fluxes on the right-hand side (RHS) 
of Eq. (1) by a new function j3 (x, t). This allows them to vary the number of 
steps in their re-initialization procedure depending on the local flow conditions, 
thus reducing artificial deformations of the interface. 

Shukla et al. in [2] assume that Eq. (1) has no physical meaning, and thus 
it can be solved in the non-conservative form without the term containing the 
interface curvature k = —V • rip. With such an assumption, the counteracting 
diffusive and compressive fluxes in Eq. (1) are projected only on the direction 
normal to the interface np 

^ = n r • V [e h \\7a\ - a (1 - a)]. (4) 

Moreover, in [2] it has been shown the key element which guarantees the suc¬ 
cessful numerical solution of Eq. (4) is the discretization of | Va| = | V-01 F (a, 7 ), 
where ip is a mapping function which smooths a (x, t) and allows to compute 
|Vet| with a smaller error. 

Since in [1, 2, 4, 5, 15, 16] the discretization and solution of Eq. (1) are 
only briefly addressed, in this paper we mainly focus on the discretization and 
solution of Eq. (1) in the framework of the second-order accurate finite volume 
method. In particular, we are interested in the case where the number of re¬ 
initialization steps in the numerical solution of Eq. (1) is N T 1 and the 
interface T is stationary. As it is described by McCaslin and Desjardins [16], 
in such circumstances T is especially prone to artificial deformations caused by 
errors in calculations of a (ipo) and np. 

The outline of this paper is as follows. In Section 2, the influence of the 
mapping function ip (a, 7 ) introduced in [ 2 ] on a convergence rate of numerical 
solutions to Eq. (1) is analyzed. In Section 3, we show that under certain condi¬ 
tions the mapping function ip (a, 7 ) approximates the signed distance function 
ipo (a) up to higher-order terms. For this reason, in the present work, it is 
proposed to use the signed distance function ipo (a) as the new mapping func¬ 
tion in the discretization of |Va| in Eq. (1). Consequently, Section 4 presents 
the selection of the mapping function ip further leads to a new, mathematically 


4 



consistent method for computation of spatial derivatives of a {ipo), and thus, 
the interface curvature. This allows us to reformulate the re-initialization and 
advection equations of the conservative level-set function a {ip o) in Section 5. 
Moreover, in Section 5, we show mathematical equivalence between the CLS 
method where the interface T is represented by a {ipo) and the standard LS 
method where the interface is represented by ipo (a), localized at the interface 
by Dirac’s delta. In Section 6 , we investigate properties of the newly formulated 
re-initialization and advection equations, in particular their convergence rates 
and errors in approximation of spatial derivatives of a {ipo) used to compute 
the interface curvature. Finally, the new re-initialization method is examined 
in several test cases with advection in order to closely inspect its conservative 
properties. 


2. Selection of the mapping function 


To assure convergence of equation (4) during integration in time r, in [2] the 
mapping function ip ( a , 7 ) that smooths a (x, t) was introduced for discretization 
of |Va|. Therein it was also noticed that ip ( 0 , 7 ) has to satisfy two conditions, 
the first condition is given by the equality 

Vq = ^ 

| Va| |V'0|’ 

as the mapping cannot change directions of the vectors normal to a (x, t) iso¬ 
surfaces. The second condition demands that the linear relation between \7ip 
and Va exists 

Va = F {a, 7 ) V-0, ( 6 ) 


where F {a, 7 ) is a known function and 0 < 7 < 1 is a constant. In this work, 
we show that the condition given by Eq. ( 6 ) can be also used to compute the 
second-order spatial derivatives of the conservative level-set function a {ip (x, t)) 


a ,ij 


F{a, 7 ) 


ip,ij + ip,iip,j 


dF~ 

da 


(7) 
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when the mapping function ip has been chosen properly, a^j = d 2 a/dxtdxj 
and i,j = 1, 2, 3. Equation (7) is required for a consistent approximation of the 
interface curvature n. 

Unlike in works [2, 4], in this paper we use the mapping functions tp for 
discretization of |Va| in equation (1). First, we note that the definition of the 
mapping function ipi (a, 7) from [2] 


Vh (a, 7) = 


(a + er 


( 8 ) 


(a + e) 7 + (1 — a + e) 7 ’ 
where originally e = 0, introduces discontinuities in the initial condition to 
Eq. (1) as it is depicted in Fig. 1(a). Two discontinuities are caused by the 
arithmetic underflow when a —> 0 and 1 — a —> 0. In order to avoid this, a 
straightforward modification of the mapping function from [2] is introduced. In 
figure 1(a) we show that setting e = 5 • 1CU 16 allows avoiding jumps in %p\. Since 
the arithmetic or floating point underflow is a purely numerical phenomenon, 
we always set e = 0 when analytical operations using Eq. (8) are performed. 

This minor modification in Eq. (8)) has a great impact on convergence of 
the numerical solution to Eq. (1). In figure 1(b) convergence of the solutions to 
Eq. (1) in the case of re-initialization of the ID regularized Heaviside function 
is presented, in this figure the distance between solutions on two different time 
levels is measured by the first-order norm 


L\ 


i 

wX 


n+1 


ai 


;=l 


(9) 


where N c is the number of control volumes and n + 1 denotes a new time level. 
In this study, the initial condition to Eq. (1) is given by Eq. (2) and we use three 
different functions in the second-order central differencing discretization of |Va|: 
a alone without smoothing, the original mapping function from [2] where e = 0, 
and the modified mapping function given by Eq. (8) where e ^ 0. As in [2, 4] 
we use the value 7 = 0.1 in Eq. (8). Re-initialization of the ID regularized 
Heaviside function is performed in the computational domain f l =< 0,1 > m 
where the interface T is located at a’r = 0.5 m; the mesh distribution is uniform 
Acc = 1 /N c , N c = 128 is the number of control volumes. At all boundaries of 
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Figure 1 : Re-initialization of the ID regularized Heaviside function H (x — 0.5): (a) 
influence of t ^ 0 in Eq. ( 8 ) on the presence of discontinuities appearing due to 
arithmetical underflow, (b) convergence of the solution to Eq. (1) with different dis¬ 
cretizations of |Vq|: a without smoothing, ipi ( 0 , 7 ) where e = 0 and ipi (a, 7 ) where 
e = 5 • 10~ 16 , 7 = 0.1 in all cases. L\ norm is defined by Eq. (9). 

the computational domain Q, the Neumann boundary condition for a (ip 0 ) is 
used. In our second-order accurate finite volume solver Fastest, the third order 
TVD Runge-Kutta method introduced in [17] is used to integrate Eq. (1) in the 
time r; the time step size is set to At = D/C 2 = eh- More details concerning 
discretization of Eq. (1) in the Fastest flow solver can be found in Appendix B. 

Since the initial condition to Eq. (1) is given by Eq. (2), we expect an 
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immediate convergence of its solution to numerical zero because equation ( 1 ) is 
initialized with its own analytical solution. However, in Fig. 1(b) it is observed 
that only the solution with e / 0 in Eq. ( 8 ) allows convergence during all 
N t = 256 time steps. 



Figure 2 : The comparison of |Va| after N T = 1 re-initialization steps of the ID reg¬ 
ularized Heaviside function with the central difference gradient approximation (c.d. 
Vq) and analytical gradient (black solid line). The |V«| in equation (1) is discretized 
using the mapping function ipi (a, 7 ) with e = 0 and e = 5 ■ 10~ 16 , 7 = 0.1. 


To explain differences in the convergence rates which are observed in Fig. 1(b), 
in figures 2 and 3 we compare the first-order derivatives of a (0o) and L\ (Va) 
norms after N T = 1 re-initialization steps. In figure 3, the L\ (Va) norm is 
defined by the equation 


Li (0) 


0 a 77. 07, 

10a 


+ e 


( 10 ) 


where 0 an , 0 num are functions calculated, respectively, analytically and numer¬ 
ically in each control volume, e = 5 • 10 -16 and 0 = Va. 

In figures 2 and 3, it is observed that both original e = 0 and modified 
e = 5 -10 ~ 16 mapping functions provide very good approximations of |Va| when 
compared with central difference gradient approximation. Differences between 
these two gradient approximations are visible only around x = 0.35 m and 












Figure 3: Distributions of Li (Va) norm defined by Eq. (10) after Ay = 1 re¬ 
initialization steps during re-initialization of the ID regularized Heaviside function. 
The error L\ is calculated for three |Va| approximations depicted in Fig. 1. 

x = 0.65 m where the jumps in i/q are present (compare results presented in 
figure 1 and in figures 2-3). If the mapping function is defined by Eq. ( 8 ) 
with e = 5 • 10 —16 , the artificial oscillations are absent since ipi is continuous 
everywhere. The continuity of ip\ guarantees convergence of the solution to 
equation (1) during all re-initialization steps as shown in Fig. 1(b). From now 
on, when referring to the mapping function ipi, we reference its definition given 
by Eq. ( 8 ) with e = 5 • 10" 16 . 

3. Relation between the mapping function and the signed distance 

function 

Figure 1(b) shows the minor modification of the mapping function ipi guar¬ 
antees convergence of the numerical solution to Eq. (1) during all re-initialization 
steps. Here, we also note that this discretization requires about N T « 140 time 
steps At to achieve the stationary solution in spite of the fact that Eq. (1) is 
initialized with its stationary solution given by Eq. (2). Therefore, the map¬ 
ping function (a, 7 = 0.1) where e = 5 • 10 -16 is not the best possible choice 
for discretization of |Vct| in Eq. (1). Hence, there is a need for an improved 
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discretization of |Va|, guaranteeing immediate convergence of the numerical 
solution to Eq. (1) towards a steady state. 

In this section, we have shown when 0 < N c 7 <1 or equivalently 0 < 7 <C 
Ax, the mapping function ipi (a, 7 ) given by Eq. ( 8 ) approximates the signed 
distance function ipo (a) up to higher-order terms. Let us first note that the 
stationary solution ( 2 ) can be rewritten as 


a W>o) 


1 

2 


1+uni (£) 


_ exp (N c i/j q) _ 

exp {N c tp 0 ) + exp (~N c ip 0 )' 


( 11 ) 


Next, we observe that 


exp (N c "ftp 0 ) 

Q/ ' - - 

[exp (N c ip, 0 ) + exp (-A^o )] 7 ’ 

_ N 7 = _ exp (—N c jip 0 ) _ 

[exp (N c ip 0 ) + exp 0 )] 7 ’ 

and we substitute Eqs. (12) and (13) into Eq. ( 8 ) to obtain 


( 12 ) 

(13) 


, = _ exp (IVc^oy) _ 

1 exp (N c i/j 07 ) + exp (-N c i[> 0 l)' 

Now, we expand exponents in Eq. (14) in the Taylor series 

A n , AT- , , ( N C.^0l) 2 , (N'cV'Ol) 3 , 

exp (N c ip 07 ) = 1 + N c i/j 0 7 H- — - 1 - - -h ... (15) 

^ A T I \ 1 AT 1 , Wd/w) 2 (IVcV’ol) 3 , 

exp {-N c ijj 07 ) = 1 - 1V c ^o7 H-^ -- b • • • (16) 

If 0 < N c 7 <C 1 is sufficiently small then the higher-order terms in Eqs. (15) 

and (16) can be neglected what gives 


exp (N c ip 07 ) = 1 + Ncipoi, (17) 

exp {-N c ij; 07 ) = 1 - iVcV’o7> (18) 


where ip' 0 
obtains 


i[>o- After substitution of Eqs. (17) and (18) into Eq. (14) one 


1 + NM 1 = 1 / J_ 
v 2 2 l 2eJ° 


(19) 


Eq. (19) is the exact relation between ipi and ^ when 0 < N c "f <C 1. Next, 
from Eq. (19) we derive the relation between the signed distance function ipo its 
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approximation ip' 0 and the mapping function %p\ (a, 7 ) 

V’o « V'o = — (2^1 - 1) • 

7 

We note that the absolute value of the gradient of ip 0 « ip'o 


dip 0 




dip'o dipi 


4 e h dip 1 

dxi 


dxi 


dip 1 dxi 


7 dxi 


( 20 ) 


( 21 ) 


since \Vipo\ = 1 is the property of the signed-distance function, see [9, 12 ]. In 
ID the second-order derivative of ipo « ip ' 0 is equal to zero 

W-(w ! ) a ^-(^)=W-(-T L )= 0 - < 22 > 

ox 1 \oxi J ox 1 \oxi J ox 1 \ 7 axi / 

The remaining question to be considered is how to select a value of the constant 
7 in the above equations; the answer to this problem is given in Section 5.1. 
Next, we show how to use the mapping functions ip 0 « ip ' 0 to compute spatial 
derivatives of a (ipo) and a (ip' 0 ), as this leads to reformulation of equation ( 1 ). 


4. Computation of spatial derivatives with mapping functions 

In this section we derive formulas for the first and second-order spatial deriva¬ 
tives of the conservative level-set function a (x, t) , see Eqs. ( 6 ) and (7). These 
new formulations exploit dependence of the level-set function a on the signed 
distance function ipo (a) or its approximation tp' 0 (ip 1 ), see Eq. (3) or Eq. (20), 
respectively. 

Let us first calculate 07 = da/dxi, i = 1, 2, 3 using Eq. ( 8 ), in this case the 
first-order spatial derivative is given by the equation 

ot,i = - ipi,i = F (23) 

7 

where C, (a, 7 ) and 6 (a) are two auxiliary functions 


p' 

II 

e 

+ 

l 

-j 

(24) 

6 (a) = a (1 — a). 

(25) 


11 



The second-order spatial derivatives of a (tpi)) are obtained directly from 
Eq. (23) and they read 


a ,ij 


S^C 2 

7 


j V’l ,ij + 

l 7 

27 ^(1 — a ) 1-7 — a 1-7 ^ 


+ (1 - 7 ) (1 - 2 a) C<5 -7 



(26) 


where a^j = d 2 a/dxidxj i,j = 1,2,3. 

Next, derivatives of a(ipo in terms of the signed distance function 

tf>o (a) are calculated, see Eq. (3), which leads to 


a 


,1 


S (a) 

e/» 


tpo,i 


(27) 


where * = 1,2,3 and (5(a) is defined by Eq. (25). The second-order derivative 
of a(t/’o ( x ,t)) is calculated directly from Eq. (27) as 


a ,ij 


S (<*) 

£h 



—(1 - 2 a) , 


(28) 


where i. j = 1, 2, 3. 

We observe that interesting similarities between Eq. (23) and Eq. (27) or Eq. 
(26) and Eq. (28) do exist when 0 < 7 <C Ax. For instance, in the case of both 
mapping functions, the right-hand sides of derived formulations are multiplied 
by 5 (a) /e^, and it is possible to show that 

1 r°° 

— / 6 (a) dipo = 1. (29) 

C/i J — 00 


Hence, S (a) /e/, approximates Dirac’s delta localized at the interface T as it has 
a compact support (on the given grid as shown in figures 2 and 3) and for e^, 7 ^ 0 
it is C°°. Consequently, 6 (a) /e^ restricts the support of a ,, a ^ derivatives to 
the region localized in the vicinity of the interface T. 

Next, we note that for 0 < 7 <C Ax (see condition required to derive Eq. 
(20)) ( (a, 7 ) « 2 inside the support of (5 (a) /eh, see Eq. (24). With this latter 
observation and from a comparison between Eq. (23) and Eq. (27) the relation 
between 7 and is obtained 


7 

4 


« e h . 
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(30) 



Eq. (30) provides the condition of equality between the first and the second order 
spatial derivatives of a (if) o) computed using the mapping function ipi (a, 7 ) or 
ip 0 (a), see also Eqs. (21) and (22). 

Finally, we recognize the relation between the present results and the result 
from the theory of distributions, see [18] page 788. When —> 0 then a (ifio) —► 
H (i/'o) and S (a) /eh —> 8 (ipo), in such case equation (23) (for 0 < 7 <C Ax) and 
equation (27) read 

VH(^ 0 ) = S^ 0 )V^ 0 , (31) 

where H (t/jo) is the exact Heaviside function and 6 (ipo) is the exact Dirac delta 
function, both functions are localized at the interface T (i /»0 ss i// 0 = 0). In par¬ 
ticular, in the direction normal to the interface rip, this definition holds only 
when | V-0o I ~ |V' 0 q| = 1 . For this reason, i/jq ~ V’o das to be the signed distance 
function to satisfy the relation between VH and S (ipo) in equation (31). 

5. Reformulation of the re-initialization equation 

After substitution of Eq. (27) into Eq. (1) the new form of the re-initialization 
equation is obtained 

^ = V - [tf(«) (|V^o| - l)n r ], (32) 

where S (a) is defined by Eq. (25) and np = V-0o/| VV’ol, see Eq. (5). We note 
that da/dr = 0 when |Vi/’o| = 1 or when S ( H ) = 0 . In the non-conservative 
form, equation (32) reads 

^ =n r -Vtf(a)(|V^|-l) 

+n r • V(|VV> 0 | - 1)8 (a) ( 33 ) 

—k (| VV’o| - 1 ) 8 (a). 

Since sign [np • VS (a) ] = —sign [V>o], the first term on the RHS of equation (33) 
resembles the term in the re-initialization equation of the signed distance func¬ 
tion introduced by Sussman et al. [10, 11]. The second RHS term in equation 
(33) contains information about the spatial distribution of a difference between 
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| V-00 (a) | and the solution to the eikonal equation |V?/>o| = 1- The third RHS 
term in equation (33) contains the interface curvature k and expresses its influ¬ 
ence on the pair of re-initialized functions a (ipo) and tpo (a). Re-initialization 
of a (i/jq) and if o (a) in equations (32)-(33) is restricted to the region of support 
of 5 (a) /eh, and thus it is localized in the vicinity of the interface T. We empha¬ 
size that the solution to equation ( 1 ) allowing re-initialization of the level-set 
function a (ipo) is mathematically equivalent to the solution to equations (32) 
or (33) allowing re-initialization of ipo {&)• 

5.1. Determination of the 7 value 

Shukla et al. in [2] and Tiwari et al. in [4] set the value of the constant 7 
in equation ( 8 ) to 7 = 0.1. They argue that this value is justified because when 
7 —> 0 the mapping function given by equation ( 8 ) tends to ipi —> 1/2 as shown 
in Appendix A. 

In this section, we investigate numerically how to select the value of the 
constant 0 < 7 <C Ax. We want to find 7 such that ipi (a, 7 ) in Eq. (20) 
accurately approximates the signed distance function ipo (a) and thus assures 
the solution to Eq. (1) with minimal error. During tests in this section, the 
number of grid nodes N c = 128 and the support width eh = Ax/2 are both kept 
constant. 

In the beginning of this study we note that Eq. (30) could be considered 
as the condition which supplies the maximal value of the constant 7 max ~ 4 eh- 
However, to derive Eq. (20), the more stringent condition 0 < 7 <C Ax is 
required; convergence of Eq. (1) with 7 = e/, is depicted in Fig. 4. 

Figure (4) presents convergence of the solutions to Eq. (1) in the case of 
re-initialization of the ID regularized Heaviside function with different values 
of the constant 7 . Errors in solutions to Eq. (1) are measured by the L 1 norm 
defined by Eq. (9). The initial condition to Eq. (1) is given by Eq. (2) which 
is its stationary, analytical solution. |Va| in Eq. (1) is discretized with the 
mapping functions ip\ (a, 7 ) and ipo (c*). In order to obtain convergence when 
7 < 10 ” 6 the size of time step At = Q t was reduced to At = eh/ 2. In the 
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Figure 4: The influence of 7 value in the mapping function ipi (a, 7 ) on convergence of 
the solution to Eq. (1) during re-initialization of the ID regularized Heaviside function. 
Li norms are defined by Eq. (9). 

present test case, the interface Xr = 0.5 m is localized exactly between the two 
neighboring control volumes xp , xf. 

When 71 = 0.1, one needs about N T « 140 time steps to achieve the sta¬ 
tionary solution. Similar to convergence with the exact signed distance function 
’ipo (cr), the most rapid convergence rate and the smallest error is obtained when 
72 = 10~ 5 . For 72 < 10 -5 , the error of the solution to equation (1) increases, 
albeit remains on a constant level. When 73 = 10 -16 , the solver needs only two 
iterations and numerical zero is achieved as shown in Fig. 4. We suspect the 
increment of the error level is a numerical effect associated with accuracy of the 
Fortran compiler (all real variables in the Fastest solver are declared in double 
precision). 

In figure 5 we compare the signed distance functions reconstructed using Eqs. 
(3) and (20) from the mapping functions ^0 (a) and ipi ( 0 , 72 ), 4> 1 (a, 73 ) after 
N t = 256 re-initialization steps. When 73 = 10~ 16 , the reconstructed signed 
distance function ip' 0 (ipi (a, 73 )) has a staircase shape due to the existence of 
numerical errors in computation of high-order roots in Eq. ( 8 ), see Fig. 5(b). 
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Figure 5: The comparison of signed distance functions reconstructed from tpi (a, 7 ) 
with Eq. (20) after N T = 256 re-initialization steps (a) 72 = 10 —5 , (b) 73 = 10 “ 16 
with tpo (a) obtained using Eq. (3). The lines with symbols are plotted every second 
or third point to improve clarity in presentation of the results. 


The latter result is empirical confirmation that the increment of Li norm levels, 
observed in Fig. 4 when 7 < 72 = 1 CP 5 , has numerical origin. In figure 5 (a), 
we observe that the signed distance function ip' 0 reconstructed from ijji (01,72) 
using equation ( 20 ) is identical with x/j 0 (o) obtained using equation ( 3 ). 
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5.2. Determination of the allowable interface width eh 

In the previous sections we have chosen th = Ax/2 to set the support width 
of 5 ( a ) /eh, as this value is also used in the literature [1, 2, 4, 5, 15, 16]. However, 
at the beginning of Section 5 it was mentioned that other than the signed 
distance function ipo (a), the Heaviside function H (ipo) is the stationary solution 
to Eq. (32) as well. Subsequently, we will demonstrate how this feature of 
Eq. (32) is preserved by the present numerical scheme which does not use flux 
limiters or higher-order essentially non-oscillatory schemes (see Appendix B for 
details). 

The lack of flux limiters in the present discretization of the re-initialization 
equation and employment of the second-order flux limiters only during the ad- 
vection step may be an advantage, as the artificial deformations of the interface 
may be avoided. A deformation of the interface due to a minmod flux lim¬ 
iter was observed in [4] during re-initialization of a (x, t) with Eq. (4). The 
interface-grid lines alignment during advection is a known deficiency in com¬ 
pressive high-resolution schemes which use down-wind to maintain sharpness of 
the interface, and switch between higher and lower order differencing schemes 
to preserve its smoothness [19, 20, 21, 22]. Moreover, the numerical artifacts 
described above cannot be accepted during the reliable implementation of phys¬ 
ical models based on the variable relaxation velocity C (x, f), and the variable 
variance eh (x, t) in Eq. (1). A good example of the physical model requiring 
abovementioned features of the numerical scheme, is the statistical model for the 
ensemble averaged description of interactions between the gas-liquid interface 
and turbulence [23, 24, 25]. 

Since in Eq. (20) we have shown that V’o ~ if' 0 for 0 < 7 <C Ax, in this 
section, for brevity, we use only i/jq for the discretization of |Va| in Eq. (1). 
In order to investigate the influence of eh < Ax/2 on convergence to the re¬ 
initialization equation, we carry out the same test as in the previous section with 
N c = 128 and with the initial condition given by Eq. (2). In the present case, the 
position of the interface is set to err = 0.6 m, and thus T is not localized exactly 
between two neighboring grid nodes xp, Xf. Additionally, eh = Ax/M where 
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Figure 6: The impact of the variable <5 (a) /th support width th = Ax/M where M < 
16, on convergence of the solution to the re-initialization equation (32). Error is 
measured by L\ norm defined by Eq. (9). 

M > 2 is an arbitrary integer number, the time step size is set to At = e^/2. 

In Fig. 6, convergence of the solution to Eq. (32) with a varying width 
of the interface e h = Ax/M can be observed. The L\ norm defined by Eq. 
(9) remains approximately constant for M = 2,4,5 and increases rapidly for 
M = 6, 8, 64. This increment, however, does not lead to the divergence of the 
present numerical solutions, we observe such behavior also when M = 16,32 
(. Li 1 < 10 -4 ). The observed increment in the L\ norm magnitude is related to 
the finite resolution of the computational grid and the selected time step size 
At. We found that errors of the solutions to Eq. (1) remain smaller than the 
truncation error (Lf 1 < 10 -16 ) when At = Ax/2 M and M = 2,..., 8. The 
explanation of this fact is straightforward if one notices that in Eq. (32) for 
eh —> 0 the value of V<5 (a) = (1 — 2a) Va increases when x ^ Xp. For practical 
reasons, later in this section we use At = eh/ 2. In such case, the errors increase 
with decreasing eh = Ax/M which follow from oscillations of the signed distance 
function ip 0 , see figures 7(b) and 8. 

In the present test case it is found that if eh = Aa:/128 then S (a) = 
a (1 — a) =0 since in the given grid \/Xi < xp : oti = 0 and Wxi > Xf : a* = 1 
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Figure 7: The influence of the variable 5(a) /th support width: th = A x/M where 
M < 16 on convergence to the solution to Eq. (32), (a) a (V’o), (b) tpo (a) after N T = 
256 re-initialization steps, the interface T is localized at xv = 0.6 m. 


where Xp,Xf are the two neighboring grid points closest to the interface po¬ 
sition xr = 0.6 m. Since 5(a) = 0 in each point of the domain only when 
a (■i/’o) = H ('i/’o) I the Heaviside function H (ipo) is the numerical solution to Eq. 
(32) as well. 

The results presented in Fig. 7 show the support width of 5(a) /th is re¬ 
stricted by the accuracy of reconstruction of the level-set functions a (i)’ o) and 
ipo (a). In order to accurately calculate gradient of a (V’o)i one needs at least four 
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Figure 8 : The impact of the variable S (a) /th support width eh = A x/M where M < 8 
on the distribution of the error \ipo/^o,an — 1|, where ipo,cm = x — xr and xr = 0.6 m, 
tpo («) is the signed distance function after N T = 256 re-initialization steps. 

points around the interface T with correctly predicted values of tpo (a). In the 
ID study presented here, this condition is satisfied when M < 4 in eh = A x/M 
(compare results in figures 7-8). 

Now, we can compare features of this new re-initialization equation with 
properties of the stationary solution to the phase-filed equation which was in¬ 
vestigated by Sun and Beckermann in [7]. The main difference lies in the fact 
that the stationary solution to the phase-field equation is always given by the 
regularized Heaviside function represented by the hyperbolic tangent profile; 
see formulation of the phase field equation in [7]. For this reason, in [7] it is 
recommend to use about 5 — 6 grid points to accurately reconstruct a (x,t). In 
the present method when eh —>• 0, the numerical solution to Eq. (32) is given 
by the Heaviside function a (ip 0 ) = H because in such case S ( H ) = 0 on 
the given grid, and thus da/dr = 0 in Eq. (32). In subsequent sections, we 
investigate how the selection of the interface width eh affects re-initialization of 
the interface T and computation of its curvature k. 
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5.3. Interpretation of results 

Let us now shortly summarize ideas introduced in the previous sections. 
Up to now the two discretizations of |Va| in Eq. (1) were introduced: |Va| = 
( 2 <5 (cc ) 1-7 /7|V0i| and |Va| = S (a) /e/j|VV’o|, see Eq. (23) and Eq. (27) respec¬ 
tively. These two discretizations are equivalent when 0 < 7 <C Ax, as shown by 
equations (20), (30) and in Fig. 5. The discretization of |Va| with Eq. (27) is 
free from round-off errors and is faster than the equivalent discretization with 
Eq. (23); in the latter case, higher-order roots of a (ipo) must be calculated, see 
Eq. ( 8 ). 

With these facts, we can now explain why equation (1) is solved with the 
largest accuracy when the mapping functions ip 0 (a) or ip' 0 (ipi (a, 72 )) are used 
for approximation of |Va|. Since ipo « ip' 0 is the signed distance function 
|V -00 (a) I ~ |V?/>q {ip\ (< 2 , 72 )) | = 1, the right-hand sides in equations (32) and 
(33) are equal to zero when the initial condition to Eq. (1) is given by Eq. (2). 
This occurs regardless of k values and explains the rapid convergence of the 
solution to Eq. (1) to a steady state as it is depicted in Figs. 4 and 6 . There¬ 
fore, the simplification of Eq. (1) to Eq. (4) put forward by Shukla et al. in 
[2] is justified only when |Va| = 6 (a) /eh\X7ipo\ where ipo is the signed distance 
function. When ipo (cc) ~ ip' Q (ipi (a, 72 )) is not the signed distance function, for 
example, due to a ( ipo ) deformations during advection, the right-hand sides in 
Eqs. (32) and (33) are not equal to zero and the re-initialization process will 
begin. 

In Section 5.2 we have shown that the stationary solution to Eq. (32) is 
also given by the Heaviside function H (ipo). This distinguishes the new re¬ 
initialization method from re-initialization of the signed distance function put 
forward by Sussman et al. [10, 11] and the solution to the phase-field equa¬ 
tion investigated by Sun and Beckermann [7]. In the ID case, this feature of 
equation (32) allows decreasing the interface width up to = Ax/4 without 
large influence on the accuracy of reconstruction of the corresponding a (ipo) 
and ipo (a) functions as shown in Figs. 7-8. The latter result suggests that for 
K dimensional problems eh = '/K Ax/4, where K = 1,2,3. In [1] and in works 


21 



that followed, the interface width is set to eh > Ax/2 and it is reported that so¬ 
lutions of the re-initialization equation (1) are not stable when eh < Ax/2. The 
new form of the re-initialization equation (1) given by equation (32) provides 
both the explanation and the solution to this problem. 

Finally, we recall when eh — > 0 then a (ipo) —► H (ipo), 8 (ck) /e^ —> 8 (ipo) 
and da/d T = 0 in equation (32). In this case, the advection equation of the 
conservative level-set function 


da (ip 0 ) , da (ipo) 


6 (a) dipo (a) 
e h _ dt 


+ Wi 


dipo (a) 
dxi 


= 0 , 


(34) 


becomes the advection equation of the phase indicator function H (ipo) which 
is discretized in the VOF family of methods, Wi in equation (34) is the i — th 
component of the interface velocity. Equation (34) shows that advection of 
the phase indicator function H ( ip 0 ) is equivalent to advection of the signed 
distance function ipo localized within the support of the Dirac’s delta function 
S (ipo)- When eh ^ 0, equation (34) describes advection of a (ipo) and advection 
of 0o (&) which is localized within the support of 8 (a) /eh- Equation (34) is valid 
in the whole domain of the solution unlike the transport equation of ipo (x, t) 
derived in [10], valid only at the interface T (wi ^ 0 when Xi € supp [<5 (V’o)] and 
Wi = 0 elsewhere). 


6. Numerical experiments 

In the following sections, we investigate the rates of convergence and as¬ 
sess numerical errors during re-initialization of both the conservative level-set 
function a (ipo) and the signed distance function ipo (a). Re-initialization is per¬ 
formed using Eq. (32) which is equivalent to Eq. (1) when |Va| is discretized 
with the signed distance function ipo or its approximation ip' Ql see Eq. (3) and 
Eq. (20), respectively. In particular, we are interested in the case when the inter¬ 
face F is stationary, = \fK/S.x/A < Ax/2 where K = 1,2,3 is the dimension 
of the problem, and the number of re-initialization steps N T ^>1. In order to 
present properties of the new re-initialization method in a broader context, sev¬ 
eral advection test cases, where equations (32) and (34) are solved alongside, are 
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performed. The numerical errors during computation of a (tfo) derivatives are 
measures of artificial deformations of the interface T due to a re-initialization 
process. In the next sections, their identification is our main concern. 

6.1. Re-initialization of stationary interfaces 

6.1.1. Regularized Heaviside function 

In figures 2 and 3 the solutions to equation (1) with the mapping function 
i/’i (a, 71 ) after N T = 1 re-initialization steps were illustrated, therein the width 
of the interface is set to = Ax/2 and the time step Ar = D /C ' 2 = e^. In 
what follows, we discuss the results obtained with the same numerical setup as 
in section (2) but after N T = 256 re-initialization steps. The results presented 
below are obtained with ipi (a, 71 ) where 71 = 0 . 1 , and with the new mapping 
functions ifo (a) and xfi ( 01 , 72 ) where 72 = 10~ 5 , see Sec. 5.1. 



x[m] 


Figure 9: The comparison of exact Va (black solid line) with its numerical approxi¬ 
mations after N T = 256 re-initialization steps of the ID regularized Heaviside func¬ 
tion. |Vck| in Eq. (1) is discretized using the mapping functions: ipo («), Vh (a, 71 ) or 
ifi ( 0 , 72 ) where 71 = 0.1, 72 = 10 -5 see Eq. (27) and Eq. (23), respectively. c.d.Va 
denotes gradient of a (x, t ) computed using the central differencing scheme. 

In figures 9 and 10, approximations to the first component of Vo (ipo) cal¬ 
culated using Eq. (23) and Eq. (27) are presented. In figure 9, it is observed 
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that both ip 0 (a) and ipi ( 01 , 72 ) reconstruct the bell-like shape of the first-order 
analytical derivative of Eq. (2) better than 1 pi ( 0 , 71 ); this result is expected in 
the light of Eq. (20). The distributions of the L\ (Vo) norms in Fig. 10 confirm 
these observations. We note that Vo is approximated with the smallest error 
in the neighborhood of the interface located at xr = 0.5 m. 



0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 


x[m] 

Figure 10 : Distributions of the L\ (Va) norms defined by Eq. (10) after N r = 256 
re-initialization steps of the ID regularized Heaviside function. Li (Va) is calcu¬ 
lated using the Va discretizations with the mapping functions: ipo (a), ipi (a, 71 ) and 
ipi (a, 72 ) where 71 = 0.1, 72 = 10 -5 , c.d. Va denotes the central differencing gradient 
a (x, t ) approximation. 

The distributions of the errors after N T = 256 re-initialization steps in the 
case of ipQ (a) and ipi (a, 72 ) are similar to the error distribution of ip 1 (a, 71 ) 
after N T = 1 re-initialization steps (compare Fig. 3 and Fig. 10). In the case 
of ip\ (a, 7 i), the L\ (Va) norm (the numerical gradient approximation within 
Eq. (10)) is varying during the integration of Eq. (1) in time r. This artificial 
deformation of Va in time r is the main reason for longer convergence of Eq. 
(1) to the steady state as shown in Fig. 4. 

Next, we discuss the accuracy of computations of the second-order spatial 
derivatives of the level-set function In Section 4, the formulas for ab¬ 

using ipi (a, 7 ) and ip 0 (a) were derived, see equations (26) and (28), respectively. 
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To compute the second-order spatial derivatives of a (ipo), one needs both the 
second and first order derivatives of the mapping functions ipo ( a ) or ipi (a, 7 ). 
In our code they are calculated using the discrete Gauss theorem equivalent to 




Figure 11 : The comparison of (a) Vr/i, (b) X7 2 i/j obtained after N T = 256 re-initialization 
steps of the ID regularized Heaviside function. 

the second-order accurate central difference gradient approximation, see [26, 27]. 

In figure 11(a) we compare the relation between Vi /71 (a) and deh/yVi/’i (ct, 7 ), 
provided by Eq. (21). One notes when t/q (a, 71 ) is used, the mapping function 
is not the signed distance function since | 4 e/ l / 7 iVi/'i (a, 71 ) | 7 ^ 1. On the other 
hand, gradient calculated using ipo (a) and i/j' 0 (ipi ( a , 72 )) gives | Vt^ol ~ |V^q| = 
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1 inside the support of <5 (a) /eh] with this, correctness of equations (20) and 
( 21 ) is confirmed. 

In figure 11(b), the second-order spatial derivatives of the mapping func¬ 
tions ipo (a) and ipi (a, 7 ) are compared. As ipi (a, 71 ) does not approximate 
the signed distance function, its second-order derivative is not equal to zero ev¬ 
erywhere in the computational domain, see Eq. (22). Since ip 0 (a) is the signed 
distance function and ip' 0 (ipi ( 01 , 72 )) is its approximation, see Eq. (20), their 
second derivatives are equal to zero when x £ supp [<5 (a) /e*,]. Consequently, 
the accuracy of approximation of the second-order derivatives using Eq. (26) 
and Eq. (28) is very similar to accuracy achieved when the first-order deriva¬ 
tives are computed with Eq. (23) and Eq. (27); compare results in Fig. 10 and 
in Fig. 13(a). 



Figure 12: The comparison of V 2 a (ip) after N T = 256 re-initialization steps of the ID 
regularized Heaviside function obtained with Eq. (26) or Eq. (28) with the second- 
order derivative calculated analytically from Eq. (2) (solid-black line) and the central 
difference (c.d V 2 a) approximation, ip is one of the mapping functions: ipo(x.,t), 
ipi (a, 71 ) or ipi ( 0 , 72 ). 

In figures 12-13, the non-zero terms in V 2 a, calculated with Eq. (26) and Eq. 
(28), are compared with analytically calculated second-order derivative of Eq. 
(2) and its central difference approximation. The first observation in Fig. 12 is 
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Figure 13: The comparison of Li (V 2 a (?/>)) norms after N T = 256 re-initialization 
steps of the ID regularized Heaviside function, ip is one of the mapping functions: 
ipi (a, 71 ), ipi ( 0 , 72 ) or ipo (a), L\ norm is defined by Eq. (10). c.d. V 2 a denotes the 
second-order derivative computed using the central differencing method, figures (a) 
and (b) present the same results but the Y-axis scales are different. 

a very good reconstruction of the second-order derivative in the case of all three 
mapping functions; see also Fig. 13 where the distribution of L\ (V 2 a) error 
defined by Eq. (10) is depicted. As expected, the second order derivatives of 
a (ipo (x,f)) calculated with \p\ ( 0 , 72 ) and ipo ( a ) are closest to each other and 
to the analytical solution. In Fig. 13, one notices when ipi (cr, 72 ) and ipo (a) 
are used, the differences between these two approximations are visible only in 
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points where the jumps of xp\ (a, 7 ) with e = 0 in Fig. 1 are present (about 
x = 0.35m and x = 0.65 m). Thus, oscillations observed in the L\ (V 2 a) norm 
can be attributed to the truncation errors due to a floating point underflow. 

In the light of the ID re-initialization studies presented above we conclude 
that the most accurate discretization of Vo) in equation (1) is achieved with 
the mapping function xjjo (a), see Eq. (27). This discretization is also the most 
natural one since the mapping between x/jq ( a ) and a (x/jq) is well defined by 
Eq. (3). We recall that after proper selection of the mapping function, Eq. 
(1) is equivalent to Eqs. (32) and (33) which are the conservative and non¬ 
conservative form of the re-initialization equation of the signed distance function 
xpo (a) and the conservative level-set function a (f/’o)- Since equation (1) is solved 
accurately when x £ supp [d (a) / eh] , there is no need for an introduction of 
additional techniques in the present solution procedure, which reconstruct the 
signed distance function xp 0 in the vicinity of the interface T . 

Influence of initial conditions and the interface width on the convergence rate. 
In most of the previous examples, the interface T was localized exactly in- 
between neighboring control volumes xp, Xf', the width of the S (a) /<q t support 
was set to eh = Ax/2 and equation (1) was initialized with its own analyti¬ 
cal solution given by equation (2). However, during advection of the level-set 
functions a (^o) and xjjQ (a), the re-initialization equation (32) must handle more 
general initial conditions. For this reason, in what follows we study influence of 
the arbitrary interface location (here Xr = 0.6m), and the support width eh on 
the convergence rate of a (xpf) and its spatial derivatives during re-initialization 
of the ID regularized Heaviside function. 

In this study, the initial condition to Eq. (1) is not given by its analytical 
solution; at r = 0 we set e/ lt o = 2 €h and we consider two widths of the interface: 
eh = Ax/2 and = Ax/4. The discretization of the ID computational domain 
is the same as described in Sec. 6.1.1. For brevity, only the results obtained with 
xpo (a) used in discretization of |Va| in Eq. (1) are presented in this section. 

In figures 14-16 it is observed that the numerical solution to Eq. (1) converges 


28 



1 



Figure 14: Convergence of the level-set function a (ipo) towards its analytical solution, 
a (x, r = 0) is set using = 2 th where (a) = Ax/2, (b) = Ax/4, the time step 

size At = e/i/2. |Vct| in Eq. (1) is discretized with ■0o (a). 

towards its analytical counterpart independent from the selected final support 
width of 5 (a) /e*,. We emphasize that re-initialization of a (ipo) with eh,o = 2 
is also possible when = Ax/M. When M > 4 the solution of Eq. (32) is 
convergent but the error level in the representation of a (i/jq) and t/>o (a) grows, 
see Fig. 17. As it is discussed in Section 5.2, this occurs due to finite spatial 
and temporal resolutions. 

Since the interface T is not localized exactly in-between neighboring control 
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Figure 15: Convergence of Va (i|io) towards its analytical solution, a (x, r = 0) is set 
using th,o = 2 eh where: (a) th = Ax/2, (b) = Ax/4, the time step size Ar = e/i/2. 

|Va| in Eq. (1) is discretized with ipo (a). 

volumes xp, xf , obtained solutions are not symmetrical as depicted in Figs. 14 
- 16. In spite of this, the shapes of the level-set function a (t/o) and its first and 
second-order spatial derivatives are correctly reconstructed. We note that in 
the limit of N T —> oo, the present numerical solution to equation (32) tends to 
its stationary analytical solution given by equation (2) and its first and second 
order spatial derivatives. The discretizations of equation (1) presented in the 
extant literature do not guarantee convergence towards its stationary analytical 
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Figure 16 : Convergence of V 2 a(i/>o) towards its analytical solution, a (x, t = 0) is set 
using th,o = 2 en where: (a) eh = Ax/2, (b) th = Ax/4, the time step size At = eh/2. 
|Va| in Eq. (1) is discretized with ipo (ab¬ 
solution, artificial deformations of the interface F that emerge when N T 1 are 
the result of this inconsistency. 

6.1.2. Circular interface 

In the case of the ID regularized Heaviside function re-initialized in Section 
6.1.1, equations (1) and (4) are equivalent since k = — V • np = 0. In 2D or 
3D cases that are discussed next, k ^ 0. Choosing the mapping function as 
the signed distance function ifo (a) allows us to write equation (1) in the form 
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Figure 17: The L\ norms defined by Eq. (9) obtained during re-initialization of the ID 
regularized Heaviside function a (ipo) with Eq. (32) and the variable interface width eh- 
The interface is located at xr = 0.6 m, at r = 0 eh, o = 2 eh, At = eh/ 2, eh = Ax/M 
and M = 2, ... , 32. 


of equations (32) or (33) and solve one of these equations without neglecting 
the influence of re on the re-initialization process. When equation (1) is solved 
correctly, equation (3) that defines the signed distance function i/jq (a) holds; 
we use this fact in the present solution procedure. 

In what follows, we assess influence of the selected mapping function and the 
interface width (e^ = Ax/2 or = \/2Ax/4) on the accuracy of approximations 
of a (ip o) derivatives, see Eqs. (23) and (27) or Eqs. (26) and (28), respectively. 
These derivatives are used in calculation of the interface curvature re according 
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to the formula 
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/1V 0o | 3 , 


( 35 ) 


which is written for the mapping function 0o (a) and is valid when Xi £ supp [5 (a) /e/,]. 
Unlike in the standard approach which uses only the signed distance function 
00, equation ( 35 ) contains terms that contribute to n exclusively away from the 
interface F. These terms are multiplied by factor (1 — 2 a) and they vanish at 
r, i.e., when a ( 0 o = 0 ) = 1 / 2 . At the interface F equation ( 35 ) reduces to k 


definition given in [12]. 

In the following sections we investigate the convergence rate of the circular 
interface curvature on five gradually refined meshes vrii = 2 4+J x 2 4+j where 
i = 1 ,..., 5 . The initial condition to Eq. ( 1 ) is given by Eq. ( 2 ) where 


00 (X, T = 0 ) 




-R, 


( 36 ) 


Li=i J 

(xo,i, £0,2) = ( 0.5 m, 0.5 m) denotes the center of the circle with the radius 
R = 0.2 m. In this test case, the computational domain is quadratic box FI =< 
0,1 > x < 0,1 > m 2 , and the number of grid nodes depends on the size of the 


grid mj. 
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Convergence of the re-initialization equation. In what follows, the convergence 
rate to the solution to Eq. (1) during N T = 256 re-initialization steps on grids 
mi, i = is presented. We compare the results obtained with two 

S (a) /eh support widths: eh = Ax/2, eh = \/2Ax/A, where the time step size is 
At = C/D 2 = eh- Unlike in the ID case, the convergence rates and L\ norms 
on gradually refined grids are practically the same when ipi ( 0 , 72 ) and ip 0 ( a ) 
mapping functions are used (compare results in figures 18(b)(d) and in figure 
4). These results again confirm the correctness of the relation given by Eq. (20). 
In the case of ipo (a), V’l (< 2 , 72 ) and for both interface widths eh the stationary 
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Figure 18: Convergence of the solution to Eq. (1) during re-initialization of the 2D 
circular interface, L\ norms defined by Eq. (9) are obtained for the mapping functions: 
ipi (a, 71 ), ipi (a, 72 ) and ipo (a), the interface width is set to eh = Ax/2 (top) and 
eh = Ax\/2/4 (bottom), At = eh, symbols correspond to every sixth or every twelfth 
iteration in time r. 

solution is achieved after about N T « 10 iterations in time r. The integration of 
Eq. (1) with the mapping function ip\ ( 0 , 71 ) requires about N T ss 25 iterations 
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to achieve the steady state solution, see figures 18(a)(c). In this latter case, the 
convergence rate is lower and the error level is higher when compared with the 
results obtained with 0 O («) and ipi (a, 72 ). 

Lack of immediate convergence of the numerical solution to a steady state (as 
in the ID case depicted in Fig. 4) is explained by additional numerical errors 
which are introduced to the solution of Eq. (1) during discretization process 
in the 2D case. The two sources of the numerical errors can be identified: 
the second-order discretization of the fluxes on the RHS of Eq. (1), and the 
computation of ip 0 (a) gradients and normals to the interface T with the central- 
differencing scheme which is known to be mathematically exact approximation 
of the spatial derivatives only in the ID case, see [26]. 

During numerical experiments it was found that differences between rates 
of convergence and their levels for the ipo («), Vh (0,72) and ipi (a, 71) mapping 
functions are more pronounced when Ax /4 < e h < Ax\[2/\ and At < D/C 2 . 
For the sake of brevity, we subsequently use the time step size At = D/C 2 = eh 
for the two different interface widths eh = Ax/2 and eh = A xy/K/4 where 
K = 2, 3 for 2D or 3D problems, respectively. 

Computation of the circular interface curvature. Next, we compute a numer¬ 
ical approximation n' of the exact curvature k using equation (35), and we 
investigate its convergence rate on five gradually refined grids after N T = 256 
re-initialization steps when eh = Ax/2 or eh = \[2Ax/A. Since a (ipo (x, t)) is 
the level-set function, k! is calculated not only at the interface cc (xr ,t) = 0.5 
but also at a(x 1 ,t) = 0.05 and at a(x 2 ,t) = 0.95, see Fig. 19. The ex¬ 
act curvature Hi of a circle on grid m; at the interface x,r is constant and 
equal to Ki = 1/7 Zi = 1 / (|x,,r — xo|). 7is the numerical approximation to 
R and is determined separately for a (V’o)> '/’o (V’l) an d 0o (ck) on each grid m^, 
i = 1,...,5, see Fig. 20. In the cases a(xl,f) and a(x 2 ,t) curvatures are 
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Figure 19: The curvature field k 3 after N r = 256 re-initialization steps of the 2D 
circular interface obtained on the grid m 3 with the mapping function ipo (a). The 
interface width is set to = Ax/2. 

defined by = 1/(7 Zi + r\) and = 1/(7 Zi + rf) where 



(37) 


(38) 


€h,i depends on the grid size and e = 5 • 10 -16 is a small constant. 

We note when the mapping function is used in discretization of Eq. (1), two 
representations of the interface E do exist. These two representations are given 
by: a (xr ,t) = 1/2 and tjj' 0 ( 1 / 2 ) = 0 or tp 0 ( 1 / 2 ) = 0 , and they are equivalent 
when eh —> 0, see Fig. 20. For this reason the LJ error is calculated as follows, 
first, the interface iso-lines are computed for each grid mf. 


1 . x“ from a(xr ,t) = 1 / 2 , 

2. xf° from if)' 0 (1/2) = 0, see Eq. (20) 

3. xf° from ip 0 (1/2) = 0, see Eq. (3). 


Later, when we refer to sets of points representing the interface iso-lines we 
use the notation x" where oj = a, i/j' 0 or i/jq. Next, the value of the curvature 
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Figure 20: Convergence of different interface representations towards the analytical 
interface (black solid line) after N r = 256 re-initialization steps. The interface width 
is set to €h = Ax/2 (top) and 6h = a/ 2 Ax/4 (bottom). The interface T is captured 
using: the regularized Heaviside function a (^o (x, t)) (solid orange dots), the signed 
distance functions 'ipQ (ipi (a, 7 )) where 71 = 0.1 (red crosses), 72 = 10 -5 (magenta 
stars), and ipo (a) (void dark blue dots) on grids: mi, m 2 , m 3 from left to right. 



Ki = 1/7 Zi = 1/max(|x^ — xo|) is computed, using equations (37) and (38) 
away from the interface T. Now, from a given field that is a numerical 
approximation of the curvature = 1/7the iso-contours given by the sets 
of points are determined on each grid ra* and for each uj. In figure 19, 
approximation of the curvature field obtained on the grid m 3 with eh = Ax/2 
is presented. 

For each grid m* the x^ and x^ ,K iso-lines are divided into Nf sample points; 
hence, the error of the interface curvature approximation is defined as 


t « _ _ 


N? 

Ik; 

1=1 


— r 


i.lh 


(39) 


where = |x“f' 


x 0 | and = |x“ z — x 0 |, x 0 is the center of the circular 
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interface. Such formula for the error accounts only for the error in the curvature 
approximation and distinguishes it from the error in the interface T position 
approximation. 

In tables 1 and 2, results of the curvature k[ convergence study on five grids 
TOj and after N T = 256 re-initialization steps are given, for the two interface 
widths eh = Ax/2 and eh = Ax\/2/&, respectively. 


m.f. 

i’l («■ 

.7i) 

V’i (« 

,72) 

V’o (x. 

,t) 

r 

a 

% (V’i) 

a 

V’o (V’i) 

a 

V’o (a) 

TOi 

4.5127e-3 

3.7965e-3 

4.0082e-3 

2.9365e-3 

4.0082e-3 

2.9365e-3 

m 2 

1.8758e-3 

6.9689e-4 

1.9332e-3 

2.1592e-4 

1.9333e-3 

2.1591e-4 

m 3 

9.6959e-4 

7.2554e-4 

8.9694e-4 

1.4100e-4 

8.9634e-4 

1.4258e-4 

7714 

8.1422e-4 

8.0377e-4 

2.6027e-4 

4.1669e-5 

2.6027e-4 

4.1669e-5 

777 5 

3.2637e-2 

3.2472e-2 

1.4827e-4 

1.0275e-5 

1.4827e-4 

1.0275e-5 

77li 

5.3185e-3 

5.2789e-3 

4.9513e-3 

4.9016e-3 

4.9513e-3 

4.9016e-3 

7772 

5.0811e-4 

5.8833e-4 

3.0678e-4 

1.5043e-4 

3.0676e-4 

1.5044e-4 

7773 

5.2778e-4 

5.8279e-4 

1.9617e-4 

1.3111e-4 

1.9617e-4 

1.3111e-4 

7774 

7.7801e-4 

7.3043e-4 

1.0271e-4 

3.7678e-5 

1.0271e-4 

3.7678e-5 

7775 

2.9551e-2 

2.9521e-2 

4.0428e-5 

1.0021e-5 

4.0428e-5 

1.0021e-5 

77li 

5.9374e-3 

6.3264e-3 

5.5758e-3 

5.8275e-3 

5.5758e-3 

5.8275e-3 

777-2 

1.7941e-3 

8.8314e-4 

1.2051e-3 

1.3091e-4 

1.2051e-3 

1.3091e-4 

777-3 

1.1891e-3 

7.3411e-4 

6.6088e-4 

1.0361e-4 

6.6330e-4 

1.0460e-4 

7774 

1.0401e-3 

7.8492e-4 

4.6755e-4 

3.6848e-5 

4.6755e-4 

3.6848e-5 

777-5 

2.6359e-2 

2.6527e-2 

2.3224e-4 

9.6314e-6 

2.3224e-4 

9.6314e-6 


Table 1 : Convergence of the curvature k'- at a (x, t) = (0.05,0.5,0.95) iso-lines (from top to 
bottom) after N T = 256 re-initialization steps, the interface width = Ax/ 2. Errors are de¬ 
fined using Lf • norm given by Eq. (39). The results are calculated on five grids mi with three 
mapping functions (m.f.) i/>i (a, 71 ), i/>i (a, 72 ), ipo(x.,t) and two interface representations 
r (a (-0o)) and T (V’o (a) ~ V’o (V’l))- 

In tables 3 and 4, the errors from tables 1-2 averaged in the narrow band 
(0.05,0.5,0.95) are presented. The results in tables 3 and 4 are depicted in 


38 



m.f. 

ip! (a. 

. 7 i) 

ip! (a 

,72) 

ipo (x. 

,t) 

r 

a 

% (V’i) 

a 

^0 (V’i) 

a 

ip 0 (a) 

m 1 

3.3388e-3 

2.2473e-3 

3.7087e-3 

4.0166e-4 

3.7087e-3 

4.0166e-4 

m 2 

1.7094e-3 

1.2542e-3 

2.5111e-3 

5.2983e-4 

2.5111e-3 

5.2979e-4 

m 3 

1.5087e-3 

1.7981e-3 

1.2191e-3 

1.6811e-4 

1.2191e-3 

1.6811e-4 

m .4 

2.3557e-3 

2.4407e-3 

3.9579e-4 

4.0151e-5 

3.9579e-4 

4.0151e-5 

m 5 

3.1094e-3 

3.0644e-3 

1.8679e-4 

1.0546e-5 

1.8679e-4 

1.0546e-5 

m. 1 

3.5091e-3 

2.6268e-3 

2.3302e-3 

1.4589e-3 

2.3302e-3 

1.4589e-3 

m 2 

8.0637e-4 

1.0710e-3 

7.4319e-4 

3.5706e-4 

7.4319e-4 

3.5708e-4 

m 3 

1.4803e-3 

1.6132e-3 

3.2301e-4 

1.4815e-4 

3.2301e-4 

1.4815e-4 

TO4 

2.8712e-3 

2.6359e-3 

1.8123e-4 

3.8763e-5 

1.8123e-4 

3.8763e-5 

m 5 

2.5595e-3 

2.5462e-3 

7.2441e-5 

1.0676e-5 

7.2441e-5 

1.0676e-5 

TOi 

6.8061e-3 

3.8563e-3 

5.4077e-3 

2.9623e-3 

5.4077e-3 

2.9623e-3 

m 2 

2.6807e-3 

1.3214e-3 

1.7051e-3 

2.7312e-4 

1.7051e-3 

2.7313e-4 

m 3 

1.8291e-3 

1.5433e-3 

7.0173e-4 

1.3536e-4 

7.0173e-4 

1.3536e-4 

TO4 

2.5678e-3 

2.3315e-3 

6.8413e-4 

3.7086e-5 

6.8413e-4 

3.7086e-5 

m 5 

3.0363e-3 

2.8712e-3 

3.1737e-4 

1.0543e-5 

3.1737e-4 

1.0543e-5 


Table 2: Convergence of the curvature at o (x, t) = (0.05,0.5,0.95) iso-lines (from top to 
bottom) after N T = 256 re-initialization steps, the interface width = Axy/1/4. Errors are 
defined using Lf i norm given by Eq. (39). The results are calculated on five grids rrii with 
three mapping functions (m.f.) ifri (a, 71 ), 'ipi ( 01 , 72 ), ipo (x,t) &nd two interface representa- 
tions T (a (V>o)) and F (ip 0 (a) « ip ' 0 (V’l))- 

Fig. 21, and they can be interpreted as the convergence rate of the circular 
interface curvature in the narrow band 0.05 < a (x, t) < 0.95 with two different 
5 (a) /eh support widths eh = Ax/2 or eh = Axy/2/\. 

In figure 21, one observes that the second order convergence rate of the curva¬ 
ture is achieved for the mapping functions ip\ (a, 72 ) and 7 > ( a ) when the signed 
distance function interface representations given by 1/0 (a) ~ ip' 0 (ip 1 (0,72)) are 
chosen. The convergence rates of the interface curvature differ dependent on 
whether the interface is represented by a(ipo) or ipQ (a). This latter observa- 
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m.f. 

V’l (a 

.7i) 

ip 1 (a. 

. 72 ) 

ipo (x 


r 

a 

V’o (V’l) 

a 

V’o(^i) 

a 

V’o 

mi 

5.2562e-3 

5.1339e-3 

4.8451e-3 

4.5552e-3 

4.8451e-3 

4.5552e-3 

to 2 

1.3927e-3 

7.2278e-4 

1.1484e-3 

1.6575e-4 

1.1484e-3 

1.6576e-4 

m 3 

8.9547e-4 

6.8082e-4 

5.8466e-4 

1.2524e-4 

5.8531e-4 

1.2613e-4 

m 4 

8.7741e-4 

7.7305e-4 

2.7684e-4 

3.8732e-5 

2.7684e-4 

3.8732e-5 

m 5 

2.9515e-2 

2.9506e-2 

1.4031e-4 

9.9756e-6 

1.4031e-4 

9.9756e-6 


Table 3: Convergence of the curvature k/ in the narrow band of 0.05 < a (xr ,t) < 0.95 the 
interface width r/ ; = Ax/2. The table contains arithmetical mean of the values in appropriate 
columns and rows of Tab. 1. 


m.f. 

V’l (« 

. 71 ) 

V’l (« 

, 72 ) 

V’o (x, t) 

r 

a 

V’o (V’l) 

a 

V’o (V’l) 

a 

V’o 

mi 

4.5513e-3 

2.9101e-3 

3.8155e-3 

1.6076e-3 

3.8155e-3 

1.6076e-3 

m 2 

1.7322e-3 

1.2156e-3 

1.6531e-3 

3.8667e-4 

1.6531e-3 

3.8666e-4 

m 3 

1.6061e-3 

1.6515e-3 

7.4794e-4 

1.5054e-4 

7.4794e-4 

1.5054e-4 

?ri4 

2.5982e-3 

2.4694e-3 

4.2038e-4 

3.8666e-5 

4.2038e-4 

3.8666e-5 

m 5 

2.9017e-3 

2.8272e-3 

1.9220e-4 

1.0588e-5 

1.9220e-4 

1.0588e-5 

Table 4 

: Convergence of the curvature in the 

narrow band 

of 0.05 < a 

(x r ,i) < 0.95 


the interface width = \f2Ax/A. The table contains arithmetical mean of the values in 
appropriate columns and rows of Tab. 2. 

tion is explained by the fact that to reconstruct the jump at F the level-set 
function a (ip o (x, t)) requires the constant number of grid points (from five to 
three dependent on selected see Fig. 7) regardless of the grid resolution 
used in simulation. At the same time, the accuracy of the representation of the 
signed distance function V’o (a), or its approximation ip q (V’l (a, 72 )), increases 
proportionally to the number of grid points N c . 

In the case of the mapping function V’l (a,ji), convergence of the interface 
curvature is not achieved although equation ( 1 ) is solved to a fully convergent 
state on all grids nii, i = 1,... ,5, see Fig. 18. This result may be explained 
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Figure 21: Convergence of the 2D circular interface curvature k' in the narrow band of 
the level-set function 0.05 < a (x, t) < 0.95 with (a) th = Ax/2 and (b) th = Ax\f2/A 
after N T = 256 re-initialization steps. In the legend a : tpo (a), the symbol on the 
left denotes the interface representation, the symbol on the right denotes the mapping 
function used to compute derivatives in Eq. (35). N c is the number of grid points in 
x, y directions. 

by the existence of a non-zero second order derivative of 'i/’i (cc, 71 ) function 
from the interface xp that impairs the calculation of the interface curvature, see 
Fig. 11(b). 

For both interface widths eh = Ax/2 or eh = \[2Ax/&, the second or¬ 
der convergence rate of is obtained for ij)a{a) and i/)\ ( 01 , 72 ), see Fig. 21. 
This confirms discussion regarding the allowable width of the S (a) /eh support 
presented in Section 5.2. The two-three grid points (see Fig. 7 and Fig. 14) 
required to reconstruct the interface curvature k with second-order accuracy 
are also needed in the construction of the flux limiters during a solution of the 
advection equation (34); the two-three grid points is a typical resolution of the 
VOF interface capturing methods [19, 20, 21, 22], 
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6.1.3. Wavy interface 

McCaslin and Desjardins [16] noticed the feedback mechanism between er¬ 
roneous normals nr and the solution of the re-initialization equation ( 1 ) when 
N t 1. This numerical phenomenon occurs due to errors introduced during 
discretization of equation ( 1 ) and leads to artificial deformations of the level-set 
function a (ipo (x, t)) increasing with time r. According to [16], this defect in the 
re-initialization procedure is particularly noticeable in regions where the inter¬ 
face T is stationary. In order to reduce this erroneous feedback between nr and 
a (if) o (x, t)) in [16] it is proposed to vary the amount of re-initialization spatially 
and to localize the solution of Eq. (1) only in the regions where the interface 
is advected or |n r • u| » 1. With this method, an additional function f3 (x,f), 
variable in time and space is introduced. The diffusive and compressive fluxes 
on the RHS of equation (1) are multiplied by /3 (x, t) to localize re-initialization 
depending on the local flow conditions. The function /3 (x, t) has to satisfy the 
condition nr • V/3 = 0 used to introduce an additional equation for /3 (x, t). We 
recognize that in the method put forward by the present paper S (a) given by 
Eq. (25) (see also Eqs. (32) and (33)), plays a role similar to the function j3 (x, t) 
in [16]. However, in the present work, we do not need to vary the number of 
re-initialization steps in time and space relative to local flow conditions. 

In this section, we carry out re-initialization of the 3D wavy interface similar 
to the test case proposed in [16]. The initial condition to Eq. (1) is given by 
Eq. (2) where the signed distance function 

i/j 0 (x, t) = 1/2 — y + A 0 sin ( 47 r;r) sin (47 rz), (40) 

defines the wavy interface, A 0 = 0.03125 m on all grids m*. Substitution of Eq. 
(40) into Eq. (2) allows us to compute the exact curvature k of the interface T. 
Convergence of k(, which is a numerical approximation of the exact interface 
curvature k,, is studied on the four gradually refined meshes: m.j = 2 4+l x 2 A+l x 
2 4+l where i = 2,... ,4 and a mesh m 3 .5 = 192 x 192 x 192 with the number 
of grid nodes in-between m 3 , m. 4 . The computational domain is a cubic box 
H =< —0.5, 0.5 > x < 0,1 > x < —0.5, 0.5 > m 3 discretized using a uniform 
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grid nodes distribution, the time step size is set to At = D/C 2 = eh- 

Using experience gained in the previous 2D test cases (see Sec. 6.1.2), two 
mapping functions if>o (a) and tpi (a, 72 ) are used for discretization of |Va| in Eq. 
(1). Similar to the previous example, we focus on convergence of the interface 
curvature in the narrow band 0.05 < a < 0.95. Additionally, as in Section 6.2, 
the influence of the S (a) jeh support width th = Ax/2 or eh = \/3Ax/A on 
convergence rate of is studied. To asses effects of the a — n r coupling on the 
interface deformations, we consider three test cases in which nr is constant or 
variable in time r: 

Ti : nr = const, Va is calculated before first iteration, N r = 1024, eh = Ax/2. 
T 2 : nr / const, Va is calculated before each iteration, N T = 1024, eh = Ax/2. 
T 3 : as the case T 2 but with = \/3Ax/4. 


Computation of numerical errors in the wavy interface curvature. In the begin¬ 
ning, we note that the curvature k derived from the analytical solution given 
by Eq. (2) and Eq. (40) is exact only at the interface T. The comparison with 
k! that is calculated away from the interface T may lead to incorrect interpre¬ 
tations of the convergence results in the narrow band of the level-set function 
0.05 < a < 0.95. This is due to the 3D wavy interface curvature Hi ^ const, in 
the direction normal and tangential to the interface T. For this reason, compu¬ 
tation of the errors in the numerical approximation of the exact curvature K t , 
is performed two ways. In the first approach, the error of the numerical solution 
is determined directly from the difference between Hi and k\, which are both 
computed in the centers of the control volumes Nf localized in the narrow band 
0.05 < a < 0.95 on the mesh rrij. The Lf, and L ^ norms on the mesh m * 
are defined by following formulas 


Nf 




1=1 


t k — 

L 2,i ~ N p 


N. 


1 1/2 


( Ki - K 'rf 


1=1 


(41) 


(42) 
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= maxO* - k[\), l = l,...,N?. (43) 

In the second approach, the convergence rate of re' is determined depen¬ 
dent on the interface representation: by the conservative level-set function 
r (cc (xj,£)) or by the signed distance function F (ip o (a)). Since in the present 
case the curvature re, is variable in space and the relation between re' and x, r 
is not known, it is very difficult to apply the procedure presented in Section 
6.1.2 for computation of re'. For this reason, we use a simplified approach using 
available post-processing tools. First, we calculate re, and re' in the centers of 
the control volumes, then the norm Lf = re,; — re' is evaluated. Afterwards, Lf is 
interpolated to the interface represented by the conservative level-set function 
T (cc (x,;, t)) or the signed distance function T (ipo (a)) obtained as iso-surfaces in 
the post-processing software. The maximal value of \Lf (T (cc)) | or |Lf (T (ipo)) \ 
on grid mt allows us to estimate L ^ i (cc) and L ^ i (ipo) norms at T (a) and 
r (ipo), respectively. 

Convergence of the re-initialization equation. In figure 22 convergence of the 
solution to Eq. (1) in the test cases Ti, T 2 and T 3 with two mapping functions 
ipo (a) and ip\ (a, 72 ) is presented. We observe that to obtain the stationary 
solution (the constant convergence level) about N T « 100 iterations are needed 
in the cases T), T 2 and about N T ss 200 iterations in the test case T 3 (here only 
ipo (a) is used). 

We note that at least second-order convergence rate of L 3 norm (defined 
by Eq. (9)) in time r may be deduced from Fig. 22 in the all three test cases. 
Surprisingly, the rate of convergence of the solution to Eq. (1) is exactly the same 
in the test cases T\ and T 2 , see Fig. 22(a)(b). Since the interface curvature re is 
implicitly included in the re-initialization equation (1) (see Eqs. (32) and (33)), 
previously discussed numerical errors introduced by calculation of gradients in 
2D-3D space have a larger impact on the convergence rate than the errors in 
computation of re. 

The convergence rate in the test case T 3 presented in Fig. 22(c) is lower than 
in previously discussed tests T), T 2 . This result is to some degree in opposition 
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n.it. N t [-] 



n.it. N t [-] 


Figure 22 : Convergence of the solution to re-initialization Eq. (1) in the test cases 
(a) Ti, (b) T 2 , (c) T 3 . |Va| is discretized using the mapping functions ipo (a) (black 
void dots) tpQ (ipi ( 0 , 72 )) (red solid dots), the results were obtained on four gradually 
refined grids mi, L 1 norm is defined by Eq. (9). 


to the convergence studies presented in figures 6 and 18. However, the wavy 
interface curvature is variable not only in the direction normal but also in the 
direction tangential to the interface T. For this reason the reduced number of 
grid points when the interface width is set to Ch = a/3Ax/ 4 < Ax/2 may lead 
to slower convergence. 

Errors in computations of the wavy interface curvature. During numerical ex¬ 
periment Ti it was found that values of the errors defined by equations (41)-(43) 
remain constant and equal to error after N T = 1 re-initialization steps, see table 
5. Similar to the convergence rates presented in Fig. 22, the values of the norms 
defined by Eqs. (41) and (43) in table 5 are identical for the mapping functions 
if 1 (a, 72 ) and x/jq (a). This result is consistent with data presented in tables 1-2 
and is expected in the light of equation ( 20 ). 
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m.f. 


ipi (cc, 72) 



V’o (a) 


L 

LI 

T K, 

^2 

T K 

^00 

LI 

T K 

T k 

■^00 

m 2 

4.6636e-2 

5.7389e-2 

1.2425e-l 

4.6636e-2 

5.7389e-2 

1.2425e-l 

ra 3 

1.1639e-2 

1.4403e-2 

3.1545e-2 

1.1639e-2 

1.4403e-2 

3.1545e-2 


5.1776e-3 

6.4156e-2 

1.4061e-2 

5.1776e-3 

6.4156e-3 

1.4061e-2 

777.4 

2.9106e-3 

3.6092e-3 

7.9169e-3 

2.9106e-3 

3.6092e-3 

7.9169e-3 


Table 5: The L£, L% and L norms obtained using Eqs. (41) and (43) in the test case T± 
after N T = 1024 re-initialization steps. 


m.f. 


ipo (a) 


L 

L\ 

T K 

T « 

^00 

7772 

4.6414e-2 

5.7396e-2 

1.2425e-l 

777 3 

1.1657e-2 

1.4478e-2 

3.1545e-2 

TO 3.5 

5.2041e-3 

6.4583e-3 

1.4060e-2 

7774 

2.9263e-3 

3.6317e-3 

7.9169e-3 


Table 6 : The Lf, L £ and norms obtained using Eqs. (41) and (43) in the test case T 3 
after N r = 1 re-initialization steps. 


In figure 23, the evolution of the errors recorded during test cases T 1: T 2 and 
T 3 are depicted. We note that after the first re-initialization step, levels of the 
all errors are equal to the values obtained in the test case T\ (compare results in 
Fig. 23 at N t = 1 with values in tables 5-6). When N T > 1 values of the errors 
grow, but after about N T ss 6 iterations they reach an almost constant level 
which remain until the end of re-initialization. Such behavior is observed in the 
three finest meshes m 2l in 3.5, 7714. The solution in the mesh m 2 is somewhat 
unresolved, since in this case we have only two grid nodes per wave amplitude 
(in [16] three nodes were used). The end values of the errors obtained in the 
test cases T 2 , T 3 are given in tables 7-8, respectively. 

In figure 24, the results given in tables 5-8 are depicted. We observe that at 
N t = 1 (and in the test case T)) errors defined by equations (41)-(43) show the 
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m.f. 


ipi (cc, 72 ) 



ipo (a) 


L 

LI 

T K, 

T K 

■^00 

LI 

T K 

T k 

■^00 

m 2 

2.7461e-l 

3.5492e-l 

1.51936 

2.7461e-l 

3.5492e-l 

1.51936 

ra 3 

1.0791e-l 

1.4944e-l 

6.2188e-l 

1.0791e-l 

1.4944e-l 

6.2188e-l 

W' 3.5 

6.8839e-2 

9.7209e-2 

3.8353e-l 

6.8839e-2 

9.7209e-2 

3.8353e-l 

7714 

5.0941e-2 

7.2103e-2 

2.7458e-l 

5.0941e-2 

7.2103e-2 

2.7458e-l 


Table 7: The LJ, and L 1 ^ norms obtained using Eqs. (41) and (43) in the test case T 2 
after N T = 1024 re-initialization steps. 


m.f. 


ipo (a) 


L 

L\ 

T K 

T « 

^00 

to 2 

2.4721e-l 

3.2123e-l 

1.210778 

m.3 

1.0624e-l 

1.4785e-l 

5.8091e-l 

TO 3.5 

6.8784e-2 

9.7189e-2 

3.7048e-l 

TO 4 

5.1173e-2 

7.2384e-2 

2.6857e-l 


Table 8 : The LJ, LJJ and norms obtained using Eqs. (41) and (43) in the test case T 3 
after N T = 1024 re-initialization steps. 


second order convergence rate. We emphasize once again that in the test case T), 
the values of all errors remain constant during all N T = 1024 steps. Hence, for 
np = const in time r the second-order convergence rate is also achieved in the 
narrow band of a (ip o)- In the test cases T 2 and X3, the first-order convergence 
rate of the interface curvature is detected in the narrow band of the conservative 
level set function 0.05 < a (ipo) < 0.95. 

As discussed previously, there is yet another way in which the errors in k' 
may be estimated. In tables 9-10, values of the L 1 ^ = max(\L K \ = |k—k/|) norms 
obtained with the mapping functions ip\ (a, 72 ), ipo (<r) in the test cases T 2 , T 3 
are given. Here, almost no difference is detected between the convergence rates 
obtained using r(a(i/)o)) or F (ip 0 (a)) interface representations, also depicted 
in figure 24(c) (f). The values in tables 9-10 show the second-order convergence 
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n.it. N x [-] 



0 200 400 600 800 1000 


n.it. N x [-] 



n.it. N x [-] 


Figure 23: Convergence of LI, L% and L norms during N T = 1024 re-initialization 
steps of the 3D wavy interface with Eq. (1) in the test cases: T 2 (left) and T 3 (right). 

| Vcr| is discretized using mapping functions ipo (cr) (black void dots) and tpb (ipi (a, 72 )) 
(red solid dots). 


rate of is detected if we take into account only the errors at the interface F (a) 
or T (compare with L 1 ^ norms presented in figures 25-27). In the case T 3 , 
the convergence rate of on the finest grid 7714 is somewhat lower; this may 
be explained by the reduced number of grid points due to the smaller interface 
width eh = y/3Ax/4: < Ax/2 used to resolve a and ip 0 ( a ). 

Finally, in figures 25-27 we can compare shapes of the reconstructed inter- 
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Figure 24: The comparison of convergence rates for LI, L% and L^ norms after N T = 1 
(void symbols) and N r = 1024 (solid symbols) re-initialization steps in the test cases 
Ti, T 2 (top) and T 3 (bottom). |Vq| in Eq. (1) is discretized using the mapping 
functions ipo (a) (dots, diamonds, squares) and ip' 0 (ip 1 ( 0 , 72 )) (triangles). 

faces obtained in the test cases T) where i = 1,2,3. The differences in the 
interface T shapes presented in these figures are barely recognizable. The 
shape of the wavy interface is preserved on all meshes used in the present study 
almost independently from the number of re-initialization steps N T . Although 
the reduction of the interface width eh leads to slower convergence of the solu¬ 
tion to the re-initialization equation, and slower convergence of its curvature n[, 
this does not affect the wavy interface shape. 

In summary, the results presented in tables 9-10 and in figure 24 confirm 
the second-order convergence rate of the interface curvature n may be achieved 
within the second-order accurate finite volume solver using the conservative 
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l k 

[1/m] 

.124 

;o.i 


.0 


I ' 01 

- 0.124 

L k 

[1/m] 

.0315 

_0.0200 


0.00 


- 0.0315 

L k 

[1/m] 

,0141 

: 0.0100 


;0.00 


--0.0100 

- 0.0141 

L k 

[1/m] 

92e-03 


-7.920-03 

Figure 25: The comparison of T (ipo) (left) and L/ = Ki — norms (right) in the test 
case Ti on four grids mi, m2, m 3.5 and mi from top to bottom. Re-initialization of 
the 3D wavy interface is performed with the mapping function tpo (a). The error Lf 
is interpolated to the interface T (ipo), the interface width is set to eh = Ax/2. 
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l k 

[1/m] 

.288 

1 


.0 


t 

-1.288 
L k 

[1/m] 

.108 

=0.08 

-0.04 

!o 

=-0.04 

j -0.08 

-0.108 

L k 

[1/m] 



-0 


"- 0.01 

I 

-0.019 

L k 

[1/m] 

9,977e-03 


I 

9 9776-03 


Figure 26: The comparison of T (ipo) (left) and Lf = Ki — norms (right) in the test 
case T 2 on four grids mi, m 2 , m 3.5 and WI 4 from top to bottom. Re-initialization of 
the 3D wavy interface is performed with the mapping function i/>o (a). The error Lf 
is interpolated to the interface T (ipo), the interface width is set to th = Ax/2 
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Figure 27: The comparison of T (ipo) (left) and Lf = Ki — norms (right) in the test 
case T 3 on four grids mi, m 2 , m 3.5 and from top to bottom. Re-initialization of 
the 3D wavy interface is performed with the mapping function i/>o (a)- The error Lf 
is interpolated to the interface T (V’o), the interface width is set to th = \/3Aa:/4. 
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m.f. 

Vh ( 

a, 72 ) 

V'o ( 

a) 

T 

a 

V’o(V’i) 

a 

ipo 

rn 2 

1.28881 

1.28783 

1.28763 

1.28772 

m 3 

1.0813e-l 

1.0771e-l 

1.0816e-l 

1.0767e-l 

TO 3.5 

1.9419e-2 

1.9067e-2 

1.9413e-2 

1.9061e-2 

7714 

1.0741e-2 

9.9776e-3 

1.0721e-2 

9.9771e-3 

Table 9: Convergence of the interface curvature 
\L^\ = \ki — k'A at the interface T (a) and T (^o) : 

measured by ' Q norm obtained 

in the test case T 2 , see figure 24(top). 


m.f. 

i>o (' 

a) 

r 

a 

ipo 

m 2 

9.3597e-l 

9.3625e-l 

m 3 

5.9071e-l 

5.9321e-2 

W 3.5 

2.0281e-2 

2.1154e-2 

7774 

1.4133e-2 

1.5113e-2 


Table 10: Convergence of the curvature re' obtained with the mapping function (m.f.) ipo (a) 
measured by /i'4. norm obtained from !/■ = re, — re' at the interface F (a) and T (V’o) in the 
test case T 3 , see figure 24(bottom). 


level-set method and the consistent re-initialization procedure. 

6.2. Tests with advection 

The primary aim of the studies performed below is verification of the new 
re-initialization method during advection of a(ipo) and ip 0 (&)• Next, we solve 
Ecp (34) with u>i = m, where it,; is the given divergence free velocity field. As we 
improve the re-initialization method first proposed in [1], similar advection test 
cases are carried out for comparison. In particular, we want to investigate the 
area (mass) conservation of the new method in the case of the variable number 
of re-initialization steps N T . Using experience gained during the tests in Section 
6.1, only the mapping function ipo (cc) is used for discretization of |Va| in Eq. 
(1) and hence Eq. (32) is solved during the re-initialization step. The interface 
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width is set to Q t = \/2Ax/A and At = D/C 2 = eh- Present investigations 
are performed in quadratic domain Q =< 0,1 > x < 0,1 > m 2 on gradually 
refined grids m.i = 2 4+ * x 2 4+i with the uniform grid nodes distribution. 

The advection equation (34) is discretized in time using the first-order im¬ 
plicit Euler method, and in space using the deferred-correction method with 
the second-order TVD MUSCL flux limiter from [28]. This type of spatial and 
temporal discretization is a default technique used in the Fastest solver for dis¬ 
cretization of advection terms in all transport equations. As the main goal of 
this paper is the improvement of the re-initialization method, the detailed solu¬ 
tions to numerical issues that arise during coupling of equations (32) and (34) 
are left for future investigations. 

6.2.1. Rotating circle 

In this section, a circular interface T revolving in the divergence free velocity 
field (ui,u 2 ) = Vo /L(y — 0.5,0.5 — x) where Vo = 1 m/s and L = 1 m is studied. 
Initially at t = 0, the center of the circle with the radius R = 0.15 m is located 
in the point (xo,yo) = (0.65 to, 0.5 to). The time step size during integration of 
equation (34) is set to At = 2n/N t , where N t = 360 • i and i = 1,..., 4 denotes 
the grid number. 

In figure 28, convergence of the two interface representations a (i/o) and 
ipo (a) towards the initial condition is presented. As previously observed in [1], 
it is clear that an introduction of the re-initialization step after the advection 
step improves the conservation of the shape and area of the advected interface. 
In the previous sections we have shown that the new re-initialization method 
does not change the position of the stationary interfaces. For this reason small 
deformations of T after one revolution may be attributed to numerical errors 
introduced by the advection scheme, see Fig. 28. The main difference between 
our results and the results presented in [1] is independence of the interface shape 
from the number of re-initialization steps N T < 32. 

To measure the convergence rate of the numerical approximation of the in- 
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Figure 28: The comparison of a (ipo) = (0.05, 0.5, 0.95) iso-lines (two columns left) and 
tpo (a) = (rj, 0, rf) iso-lines (see Eqs. (37) and (38), two columns right) obtained after 
one revolution of the circular interface on grids rm, i = 1,. .., 4 (from top to bottom) 
with the initial condition (black solid line). The results were obtained with N T = 0 
(orange double-dashed line), N T = 1 (magenta dashed-dotted line) and N T = 32 (dark 
blue dashed line) re-initialization steps. 


terface T towards analytical solution, the first order norm is introduced 

Nf 


t r _ _ \ ' | num „ext 

L l,i - T7J 2^ l r i,7 r 


Li I) 


(44) 


i 


where Nf is the number of points in the interface T representation given by the 
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iso-lines of a (^o) = 1/2 or (a) = 0, r”“ m = |x™ m - x 0 |, iff = |xff - x 0 |, 
where x"“ m and xff are points obtained from the numerical and the exact 
approximation of the interface T (given by the initial condition) on the grid 
rn t . As in the present numerical scheme, we use the first-order discretization 
method in time. This convergence rate of T towards the initial condition is 
also observed in Tab. 11 for both interface representations. The convergence 


T 


a (ipo) 



ipo (a) 


N t 

i 

16 

32 

1 

16 

32 

mi 

8.3833e-3 

8.1417e-3 

7.9349e-3 

8.5084e-3 

8.3689e-3 

8.1521e-3 

m 2 

4.4U9e-3 

4.5067e-3 

4.4682e-3 

4.3593e-3 

4.4516e-3 

4.4166e-3 

ra 3 

2.1589e-3 

2.2262e-3 

2.2192e-3 

2.1621e-3 

2.2296e-3 

2.2221e-3 

7714 

1.0916e-3 

1.1243e-3 

1.1239e-3 

1.0923e-3 

1.1247e-3 

1.1244e-3 


Table 11: The first-order convergence of LJ ^ norm, computed after one revolution of the 
circular interface represented by T (a (V’o)) and F (ifto (a)) with N T < 32 re-initialization steps. 


rate of the reconstructed interface T towards the initial condition is related to 
details in the coupling between equations (32) and (34) and could be improved 
by introduction of the liigher-order discretization schemes. We note that the 
values of L\ i norm in table 11, remain almost the same in spite of a different 
number of re-initialization steps N T used during simulation. 

Since the present method uses two representations of the interface a (ip o) and 
ipo (a), the area (mass) conservation can be estimated in two different ways. In 
order to distinguish between these two possibilities, a cumulative error is defined 
as 


N t 


Et = — ^2 E n where 


En = \S n num - S n ext | = J \a n (iPo) - aT I dS, 


(45) 


and a ext , ipQ Xt are defined by the equations (2) and (36), respectively. N t in 
equation (45) denotes the total number of time steps on the grid m, during one 
revolution of the interface T, E n is a difference between the numerical approxi- 
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mation of the surface S™ um and the exact surface S r p xt which are calculated on 
each time level n. The instantaneous position of the circular interface center 
(x c i Vc) required to compute of S r e l xt is obtained from 


x c 


^ Lii,j _ X) i,j a i,jyi,j 

> Vc 

i,j / jjj &i,j 


(46) 


We found that convergence rates of the area (mass) depend on the region 
of integration of a (?/’o) hr Eq. (45). If the region r\ = {Xi | 1 — a ( ipo ) > 0.5} 



0 n/2 it 3n/2 2 n 

t[s] 

Figure 29: The area conservation calculated in the regions: (a) ri = 

{xi\l - a (ipo) > 0.5}, (b) r 2 = { Xi | ipo (o) < 86^} during one revolution of the cir¬ 
cular interface with N T = 1 or N T = 32, the error is normalized with S an = tiR 2 . 
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Figure 30: Convergence of the cumulative error Et defined by Eq. (45) with the different 
number of re-initialization steps N r <32. E t is calculated in the regions: (a) r\ = 
{xi\l - a (tp 0 ) > 0.5}, (b) r 2 = { x / | i/’o (o) < 8r/; } and is normalized with San — 
■kR 2 . The boxes on the right, present convergence during last ten time steps At, the 
number and position of symbols is arbitrary. 


or ri = {cci | ipo (a) < 0} is chosen, then the first-order convergence rate of area 
is obtained, see figures 29(a)-30(a). Similar convergence rates of the Heaviside 
function were reported in [1] although therein the second-order discretization 
in time was employed. When the integration is carried out in the region r 2 = 
{Xi (a) < 8e/j} or in the whole computational domain Q, then almost the 
second-order area convergence rate may be deduced from the results depicted in 
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figure (29) (b). In this latter case, convergence of the interface at each time step 
At with N t < 32 (see figure 30(b)) is less evident but may be also observed for 
example on the grid 7714 . The differences in the magnitude of E t errors visible in 
figures 30(a) (b) may be explained by the more accurate interface representation 
with the signed distance function, and thus smaller error E n during the area 
(mass) computation in the region ?’ 2 , see Eq. (45). 

Here and in the following examples, it becomes clear that S num calculated 
in the region ri = {a;* | 1 — a (ipo) > 0.5} indicates how the area (mass) varies 
during topological changes of the interface (stretching, break up, coalescence), 
whereas S num calculated in the region r2 = {#* | 1/0 i a ) A 8 e/, } allows examina¬ 
tion of the total area (mass) conservation. Since the circular interface is rotated 
as a rigid body in the present case, the errors integrated in both regions ri and 
72 remain on constant levels, see Fig. 29. 

During rotation of the circular interface without deformations, obtained 
area convergence rates are almost independent from the number of selected 
re-initialization steps N T (compare results with different N T in Fig. 29 and in 
Tab. 11). The results presented in figure 30 suggest that in most cases up to 
four re-initialization steps should be sufficient to preserve the shape and area of 
the interface T. 

6.2.2. Vortex test 

To test the new re-initialization method in a more complex velocity field, we 
use 


Mi = —Vo sin 2 ( kx ) sin ( 2ky ), U 2 = Vgsin 2 ( ky) sin ( 2kx ), (47) 

where k = n/L, L = lm and Vo = 1 m/s (similar to the test case in [1]). A circle 
with the radius R = 0.15 m is initially located at (x 0 , yo) = (0.5 m, 0.75 m). The 
simulation time is t = 2 s, as after time T = Is the flow field is reversed so the 
exact solution should be obtained after the same number of time steps. In order 
to satisfy the Courant condition Cr ss 0.2, the time step size is set to At = Ax/8 
to give the total number of time steps per revolution N t = 8T/Ax = 8TN C . 
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The iso-lines of ipo ( a ) obtained on grids to*, i = 1,... ,4 with N T = 2 or 
N t = 16 are presented in figure 31. For the sake of brevity, iso-lines of a (ipo) are 
omitted as they are almost identical with iso-lines of i/jq (a). The most important 
observation in Fig. 31 is that on finer grids, the impact of the number of re¬ 
initialization steps N t on the interface shape becomes negligible. This result 
shows convergence of the new re-initialization method and is in agreement with 
studies presented in Fig. 28. 

When the present results are compared with the results presented in [1], one 
notices due to the smaller interface width in our simulations eh = \f2Ax/k < 
Ax/2 our results are qualitatively similar to the results from [1] on a grid twice 
as fine. This illustrates the interface width eh is a very important parameter in 
the present method; eh decides not only whether the interface curvature may 
be calculated (see Section 5.2) but also governs the topological changes of the 
interface when the grid resolution is not sufficient to resolve it, see Fig. 31. 


r 


a('0o) 



ipo («) 


N T 

2 

8 

16 

2 

8 

16 

mi 

2.1424e-l 

2.3195e-l 

6.8181e-l 

2.1434e-l 

2.3175e-l 

6.8196e-l 

m 2 

7.3558e-3 

1.0785e-2 

1.1301e-2 

7.3921e-3 

1.0773e-2 

1.1305e-2 

m 3 

3.2984e-3 

3.4034e-3 

3.5041e-3 

3.3082e-3 

3.3917e-3 

3.4896e-3 

ra4 

1.6404e-3 

1.6512e-3 

1.6444e-3 

1.6467e-3 

1.6542e-3 

1.6433e-3 


Table 12: Convergence of Lj i norm, computed at the end of the vortex test with the N T < 16 
re-initialization steps. The circular interface is represented by T (a (ipo)) and T (ipo (a)) iso¬ 
lines, compare with results in figure (32). 

After rotation of the interface in the opposite direction with the same num¬ 
ber of time steps, the circular interface F is reconstructed with the first order 
accuracy, (see results in figure 32 and in table 12). Similar to the case without 
interface deformation, the errors in Tab. 12 do not change much with N T . Con¬ 
vergence towards the initial condition may be observed on the gradually refined 
grids rrii, i = 1,..., 4. 
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Figure 31: Vortex test on four different grids rrn, i = 1,..., 4 from left to right at times 
t = 0, t = 0.5, t = 1.0, t — 2.0 from top to bottom. Figure depicts ipo ( a ) = (rf, 0, rf) 
iso-lines (see Eqs. (37) and (38)) with N T = 2 (magenta, dashed-dotted line) and 
N t = 16 (dark blue, dashed line) re-initialization steps. 


The L\ i norm defined in Eq. (44) is not the best measure of interface depar¬ 
ture from its original shape as it does not detect oscillations in the reconstructed 
interface. Closer investigations of the results obtained on the grid rri 4 (dark blue 
dashed line in Fig. 32) reveal that the regularity of the final interface shape at 
t = 2s grows with the N T number. Since re-initialization improves ipo (a), the 
errors of rip remain smaller when N T is larger and for this reason oscillations 
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Figure 32: Convergence of tpo (a) = 0 iso-lines on grids rm, i — 1,..., 4 towards initial 
condition (black solid line) in the vortex test case with (a) N r = 2 , (b) N T = 8, (c) 
N t = 16 re-initialization steps. 

does not appear in the interface shape reconstructed at the end of simulations 
(compare results on the grid 7714 in figures 32(a) and 32(b)(c)). 

In figure 33, the area conservation in the vortex test case is presented. Like 
in rotation of a rigid body, we also calculate (integrate) errors in two regions ry 
and r 2 in the present case. Since in this case the circular interface is strongly 
deformed through the vortex velocity field, large variation of the error integrated 
in the region ry is observed, see Fig. 33(a). As the area of the interface F 
increases when it is deformed, the value of the error becomes smaller and then 
returns to its initial level at the end of the simulation. This result confirms the 
present re-initialization method is conservative. Similar variation may also be 
observed in the case of the error obtained in the region r 2 , however the effect is 
very small, and therefore cannot be observed in Fig. 33(b). The distributions of 
the errors in Fig. 33(a)(b) confirm that at least the first-order convergence rate 
of the area (mass) is achieved by the present numerical method. Moreover, we 
note that the number of re-initialization steps N r performed on each time step 
At has a small impact on the calculated area when the mesh is sufficiently fine 
(compare results with different N T in figures 31 and 33). 

To investigate behavior of the new numerical method during longer inte¬ 
gration times (here t = 3 s), we perform the vortex test without reversing the 
velocity field after the first revolution. In this test case, we stretch the inter- 
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Figure 33: The area conservation calculated in the regions: (a) n = 

{ay | 1 — a (tpo) > 0.5}, (b) r 2 = {xi | ipo (a) < 8th } during vortex test with N T < 16. 
Error S„ urn defined in Eq. (45) is normalized with S an = kR ?, the number and position 
of symbols is arbitrary. 


face until it is broken into bubbles due to insufficient number of grid points 
required for its reconstruction. As our advection and re-initialization methods 
are conservative, and the re-initialization method keeps the prescribed thickness 
of the interface constant, this is the only possible scenario which may occur 
when the grid resolution is insufficient to resolve the interface T. We emphasize 
that there is no physical mechanism behind the abovementioned break up of 
the interface. This feature of the present method is direct consequence of its 
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Figure 34: ipo (a) = 0 iso-lines, in the vortex test on grids rrn, i = 2,... ,5 (from top to 
bottom) after T = 0, T = 1, T = 2, T = 3 revolutions (from left to right), the width 
of the interface th = \/2Ax/4, N T = 4 on each time step At, At = eh- 

conservative properties. 

Using experience gathered in the previous examples, we use N T = 4 re¬ 
initialization steps in each single time step At, e/,. = \/2Ax/4 and Ar = eh- 
The results from this study, obtained on four grids mi where i = 2,... ,5, at 
four different time moments T, are presented in Fig. 34. In order to visualize the 
interface thickness on each grid m*, at T = 0 three iso-lines of ipo (ct) = {rj, 0, rf) 
(see Eqs. (37) and (38)) are depicted. 
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Figure 35: The errors in the area conservation calculated in the regions (a) n = 
{xi | 1 — a (tpo) > 0.5}, (b) v 2 = {xi | ipo (a) < 8 eh} during vortex test with N T = 4 
on each time step At. Error Snum defined in Eq. (45) is normalized with S an = nR 2 . 

We note that the longer the integration time is, the thiner the interface F 
filaments become. However, if the mesh resolution is insufficient to resolve the 
interface F, then the interface break up will occur as this is the only way to 
preserve the area (mass) correctly. 

Errors in the area conservation are depicted in figure 35, where, as in the 
previous examples, S num is computed by integration in two regions t*i and r 2 - 
The results in Fig. 35 show that the area is conserved in spite of large defor- 
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mation of the interface T on all grids m*, i = 2,..., 5. As discussed previously, 
the error calculated in the region r i permits tracing the impact of the interface 
deformation on the area (mass) conservation. We note that after break up of 
the interface into separate bubbles on grid m 2 at the time about T = 2.5 s, val¬ 
ues of S num integrated in region r\ cease to drop as the interface deformation 
is stopped (see Fig. 35(a)). It is expected after large enough integration times 
t 3 s, at each grid, the interface will finally disintegrate into droplets that 
can can be transported in the given velocity held. 

The total error in the area (mass) conservation, obtained by integration of 
a (if) 0 ) in the region r 2 , is almost constant on all meshes used in the present 
study. The loss of total area (mass) is present but negligibly small as shown in 
figure 35(b). 

7. Conclusions 

In this work, a new re-initialization method of the conservative level-set func¬ 
tion was introduced and verified. We have shown that the re-initialization and 
advection equations of the conservative level-set function a (ipo) are mathemat¬ 
ically equivalent to the re-initialization and advection equations of the localized 
signed distance function ip 0 (a) (see equations (32) and (34)). It was also proven 
that the RHS of re-initialization equation (32) is equal to zero when |V0o| = 1) 
and a (ipo) — H (ipo)- These two solutions to the re-initialization equation 
(32) permit computing the interface curvature, and controlling its thickness 
eh < Ax/2 which assures the improvement of the spatial resolution of the new 
method, see discussion in Sec. 5.2 and the results in Sec. 6. 

In the present re-initialization procedure, the conservative level-set func¬ 
tion a (ip 0 ) remains continuous and bounded inside the support of 6 (a) /e^, = 
a (1 — a) /eh function. For this reason, we do not need to use the fast marching 
method, required in other re-initialization techniques reported in the literature, 
in order to extend the signed distance function ip 0 (a) away from the interface 
r. The continuity of the solution assured by the new re-initialization method 
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avoids artificial interface deformations, which are eventually introduced by the 
flux limiters typically used to bound counteracting compressive and diffusive 
fluxes on the RHS of Eq. (1). In the light of the other results presented in the 
literature, the new re-initialization method shows fast convergence and in the 
ID cases, reconstructs exact behavior of the analytical solution to the partial 
differential equation (1), see Fig. 4. 

Discretization of the first and second-order spatial derivatives of a (ip o), 
which is consistent with the discretization of the re-initialization equation and 
the theory of distributions, achieves the second-order convergence rate of the 
interface curvature k, see Eq. (35). Such level of accuracy is obtained in the 
second-order accurate finite volume solver without geometrical reconstruction 
of the interface or an introduction of the higher-order essentially non-oscillatory 
interpolation schemes in the discretization procedure. 

The advection tests performed in Section 6.2 show the new method conserves 
the total area (mass) with almost second order accuracy, and the shape of the 
advected interface is independent from the number of the re-initialization steps 
if the spatial and temporal resolutions are sufficient. 

The new re-initialization method presented herein may be considered a back¬ 
bone of the conservative level-set method, a potentially good replacement for the 
compressive interface capturing schemes commonly used in the fast multiphase 
flow solvers. Its implementation is simple and its strengths include accurate 
reconstruction of the interface curvature and conservation of the area (mass) 
without smearing the interface. 
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Appendix A. Limits of the mapping function 


In this appendix, the limits of the level-set function a(ipo (x, t)) and the 
mapping function ip i (a, 7 ) given by Eq. (2) and Eq. ( 8 ) are calculated. Since in 
this section we consider analytical properties of ipi (a, 7 ) next e = 0 in equation 
( 8 ). We assume that the number of grid nodes N c —► 00 and 0 < 7 <C Ax is a 
small constant. The analytical solution given by Eq. (2) may be rewritten as 


1 


1 + tanh 


2e?i 


1 

1 + exp (—2N c ip 0 ) 


(A.l) 


where N c = 1/Ax and ipo is the signed distance function given by Eq. (3). 
Additionally, we rearrange Eq. (14) into 

exp {N c ip 07 ) 1 


i>i = 


(A.2) 


exp {N c ip 0 j) + exp (-jV c 0 O 7) 1 + exp (-21V c 0o7)' 

First, we note that if 7 —>• 0 in Eq. (A.2), ipi —> 1/2 as it was observed in [2]. 
Let us now consider three limits of Eq. (A.l) and Eq. (A.2) when N c 7 —► 00 
and 0 < 7 < Ax: 


1 . when 00 < 0 

2 . when 00 = 0 

3. when 0o > 0 


ip! — y 0 , (x — y 0 , 
ip! = 1/2, a = 1/2, 
0i —> 1, a —> 1 


Hence, when the number of grid points N c —> 00 and 7 = const functions 0i and 
a become the Heaviside function H (ipo). The phase indicator function H(ipo) 
is typically discretized in the volume of fluid (VOF) methods [ 8 ]. 


Appendix B. Discretization of the re-initialization equation 


In this appendix, the discretization of Eq. (1) or equivalent Eqs. (32) and 
(33) in the framework of the second order accurate finite volume method is 
presented. Since all terms on the RHS of Eq. (33) are in the form of the mass 


forces, after its integration in the control volume Vp one would obtain 

g| p = [n r .Vd(a)(|V0|-l)]|p 

+ [n r -V(|V0| —l)<5(a)]| P (B- 1 ) 

- [k(|V0| - 1 )6(a)] | P , 
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where nr = V</>/|V^>| and <j> is equal to a (x, t) or is one of the mapping functions 
ipi (a, 7 ), ipo («)■ Eq. (B.l) is the discretization to the non-conservative form of 
Eq. (1), all terms with sub-script P can be obtained from the values stored in 
the centers of the control volumes. 

In the present work, we use the conservative discretization to Eq. (1) or 
equivalent Eq. (32). After integration of Eq. (32) in the control volume Vp, 
employment of the Gauss theorem and mid-point rule at the centers of faces / 
and in the centers of control volume P one obtains 


da 

dr 


p 


1 nb 

y~ P ( a ) (I V '^ ) l - n r ' n ]/ S f , 
P /=1 


(B.2) 


where nt, is the number of neighbors of the control volume P, [nr • n]^ is dot 
product of the normal nr = V/>/|V(/| interpolated to the face / and normal n/ 
at the face /, <5 (a) is defined by Eq. (25); |V</| is computed using the second- 
order central-difference approximation to the <f> gradient components, at face 
/ = e this approximation reads 


dcj) 

dx 1 e 

d(f> 

dx 2 e 

dcj) 

dx 3 e 


{<f>E ~ <t>p) 

Ax 

(< i>N + (j>NE ~ (j>S — <t>SE ) 
4 Ay 

{(j>T + 4>TE — <j>B — 4 >Be) 

1A: 


(B.3) 


where subscript E, N,T,... represent the centers of the neighbor control vol¬ 
umes. The interpolation to the faces f of the control volume P is performed 
using second-order accurate linear interpolation scheme, which on uniform grids 
simplifies to </>/ = 1/2 {(f>p + <f>p), see [26, 27]. Dependent on the test case, <j> in 
Eq. (B.3) is calculated using a (x, t) or ij) obtained from the mapping functions 
given by Eq. (3) or Eq. ( 8 ) with different 7 values. 
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